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ABSTRACT 

This investigation is concerned with minimum weight designs of arch 
structures. Using the finite element method, the arch is modeled by 
contiguous bar-beam elements. Element stiffness coefficients in terms of 
local degrees of freedom are related to system degrees of freedom through 
local to global coordinate transformations. After coordinate transformations, 
element stiffness coefficients are assembled into FEM stiffness equations for ’ 
the arch structure. An objective function for weight minimization, with 
constraints on failure, arch geometry, and section dimensions, is minimized 


by the DOT optimization code. Results are presented for a number of cases. 


iv 





TABLE OF CONTENTS 


TC INTRODUCTION: wiccni us teioktasehn eee marcy see ae eot eames 1 
Il. PROBLEM FORMULATION .......... cc cee eee cece eee e eee eee 8 
A. THE OBJECTIVE FUNCTION ........... 0... cece cece cece eee 9 
B. THE STRENGTH CONSTRAINT .............. 0. cee ee eee ees 10 
C. GEOMETRIC CONSTRAINTS ............ cece cece cee ee eeee 11 
D; “SIDE CONSTRAINT S © sos ont 5dies woieedd ween tadiu eee 12 
E. OPTIMIZATION SOFTWARE ........... cc cece cece ee eee eee 12 
TE STRESS: ANALYSIS 65s anes Hide avec teks weosaseiineene ts ink 15 
IV. FINITE ELEMENT ANALYSIS ............... 22 ec cece cee eeee 18 
A. THE BEAM EQUATION DEVELOPMENT ................666: 18 
B. THE BAR EQUATION DEVELOPMENT ................eee0e: 21 
C. THE ELEMENTAL STIFFNESS MATRIX ................2000- 23 
D. COORDINATE TRANSFORMATION ............. 2c ce ee eecees 27 
E. THE ELEMENTAL SYSTEM OF EQUATIONS ................- 30 
F. THE GLOBAL SYSTEM OF EQUATIONS ...............-..00- 30 
V. PROGRAM DESCRIPTION AND CAPABILITIES .................. 33 
VI. PROGRAM VERIFICATION ............ 2c cece cece cece ccc eeeee 39 
V1 ‘CASE STUDIES ies cowie dad cow seeasiewi estan bese as 47 
A. CASE 1: THE OPTIMIZED FRAME ................. eee eees 48 
v 


B. CASE 2: THE OPTIMIZED CANTILEVER ARCH .............. 50 
C. CASE 3: THE OPTIMIZED CANTILEVER ARCH WITH AXIAL END 


DOBD cipro Ss stat Se Sea wean pi GSE PR DASS EWAE aa He Rees 51 
D. CASE 4: THE CANTILEVER ARCH UNDER A DISTRIBUTED 
LOAD (os cen cakes eases setae ae Rees ee bee wR we eaatey 54 
E. CASE 5: THE SIMPLY SUPPORTED ARCH ................... 56 
F. CASE 6: THE FIXED-FIXED ARCH .............. ccc scenes 58 
VHT CONCLUSIONS ooicadeea owas teens hp SESS Saw ate been eaaees 61 


APPENDIX A JUSTIFICATION FOR OMITTING SHEAR STRESSES .. 63 


APPENDIX B DOT USERS MANUEL, SECTION 2.1 ................. 66 
APPENDIX C DOT USERS MANUEL, SECTION 2.2 ................. 67 
APPENDIX D VERIFICATION AND CASE STUDY OUTPUT .......... 71 
APPENDIX E ARCH_OPT.FOR SOURCE CODE ................4... 92 
LIST OF REFERENCES 54. ca:s icine eee esein wesw eau eee 107 
INITIAL DISTRIBUTION LIST ........ ec cece cece cece cece enee 108 





Tables 


TABLE 5.1: ARCH_IN.DAT FIELD PARAMETERS ................... 35 
TABLE 6.1: VERIFICATION PROBLEM #1 SUMMARY OF 
DISPLACEMENTS. :3 cis csinsagae' ve eiwiceus yeatacwves wes 40 


TABLE 6.2: VERIFICATION PROBLEM #1 SUMMARY OF STRESSES .. 41 
TABLE 6.3: VERIFICATION PROBLEM #2 SUMMARY OF RESULTS ... 43 
TABLE 6.4: VERIFICATION PROBLEM #3 SUMMARY OF 


DISPLACEMENTS: o:s.0saissca tes sce oses sesso rectewaacine 44 
TABLE 6.5: VERIFICATION PROBLEM #3 SUMMARY OF STRESSES .. 45 
TABLE 7.1: CASE 1 SUMMARY OF RESULTS ................000005 50 
TABLE 7.2: CASE 2 SUMMARY OF RESULTS ................00000- 51 
TABLE 7.3: CASE 3 SUMMARY OF RESULTS ...............-..005. 54 
TABLE 7.4: CASE 4 SUMMARY OF RESULTS ................-..05- 56 
TABLE 7.5: CASE 5 SUMMARY OF RESULTS ...................--. 58 
TABLE 7.6: CASE 6 SUMMARY OF RESULTS ...............-2..06- 60 

vii 





Figures 


Figure 1.1: Pont du Gard of the Nimes Aqueduct .................-2000- 3 
Figure 2.1: Bar-Beam Model of the Arch ........... 000 ccc ecececceeees 9 
Figure 2.2: Arch Dimensions ............ 0c ccc cece cece cece ees eecens 10 
Figure 3.1: Normal Stresses due to Bending and Axial Force ............ 16 
Figure 4.1: Beam and Bar Element Degrees of Freedom ............. «. 24 
Figure 4.2: Beam-Bar Element Degrees of Freedom ................--. 25 
Figure 4.3: Displacement & Force Coordinate Transformations .......... 28 
Figure 5.1: ARCH_OPT Program Structure ............ 0.2 ceeeecsees 34 
Figure 5.2: ARCH_IN.DAT Variable Implementation .................. 38 
Figure 6.1: Verification Problem #1 ............ 0. cece cece eee ceceeees 40 
Figure 6.2: Verification Problem #2 ............ cc ccc ec cece cece cecees 42 
Figure 6.3: Verification Problem #3 ............. cece cece cece eee e ewes 44 
Figure 6.4: Arch Verification Problem #4 ..............cccecccscccees 45 
Figure 7.1: Methods of Transferring a Load ..............00eeeeceeees 48 
Figure 7.2: Case 3 Problem Geometry ............. 2. cece cee ceseeees 52 
Figure 7.3: Moment Diagrams of Case 2 and Case 3 ...............005- 53 
Figure 7.4: Case 4 Problem Geometry ............... cece cece eee eeee 54 
Figure 7.5: Case 5 Problem Geometry .............. 00. e sce e ence neces 57 
Figure 7.6: Case 6 Problem Geometry ............. 0... cece eee eeeeees 59 





D 
DOT 


Table of Symbols 
the cross-sectional area 


width of the i element 


pari pietance from the center line to the outermost fiber of the element; 
c= 


the Domain of the problem 
Design Optimization. Tool software from VMA Engineering 


Young’s Modulus of the arch material 
the bar-beam force vector in the global coordinate system 


the bar-beam force vector in the elemental coordinate system 

the bar elemental force vector 

the beam elemental force vector 

Concentrated axial force 

the bar system force vector 

the beam system force vector 

the bar system force vector including the boundary term vector U 


the beam system force vector including the boundary term vectors M 
and VY 


Finite Element Method 
the column vector of linear shape functions 


height of the i element 


cross-sectional moment of inertia 
the bar-beam elemental! stiffness matrix in x-y coordinates 


the bar-beam elemental stiffness matrix in local coordinates 








a a 


= 


28 


“4 


Pe Me y 


the bar elemental stiffness matrix in local coordinates 
the beam elemental stiffness matrix in local coordinates 
the bar-beam system (global) stiffness matrix 

the bar system (global) stiffness matrix 

the beam system (global) stiffness matrix 

length of the i** element 

the total length of the given structure 

the differential operator 

Moment 

Maximum Moment 

Concentrated Moment 

the moment boundary term vector 

the total number of elements 

axial loading 

lateral loading 

concentrated load 

the column vector of cubic shape functions 


the ratio of the maximum shear stress to the normal stress due to 
bending; r = t,_,/s, 


the radius of the arch 

the Residual function 

the ceater-line coordinate of the arch 
yield strength of the arch material 


axial displacement 
the approximate axial displacement 





l=] 


<< 4Q<¢«iea 


a 8 PF BoD HM Yar ROM & 
a 


F 





the vector of axial displacements 
the bar equation boundary term vector 
lateral “displacement” 


the approximate lateral “displacement” 
the vector of lateral displacements and siopes 


the shear force 


the snear force boundary term vector 


the horizontal axis 


the vertical axis 


the zero vector 
the angle the i** element makes with the x-axis 
the perpendicular compliment of «, 


the bar-beam displacement vector in the x-y coordinates 
the (6x1) bar-beam displacement vector associated with k’’ 


the exact analytical solv tion 

the (6x6) local transformation matrix 

the subtended arc of the arch 

the normal stress due to bar (axial) behavior 


the normal stress due to beam (bending) stress 
the maximum stress developed in the it element 
the total normal stress 


the maximum shear stress 





Acknowledgment 


My thanks and gratitude to Professor Dong Soo Kim for his time and 
effort in teaching me the ways of DOT and to Professor Salinas for his open 
door and infinite patience. 


I, INTRODUCTION 


The arch has been employed in the fabrication of engineering 
structures for over five thousand years. Its suitability to compressive load 
design made it especially favored when masonry, not steel, was the principle 
building material. (Due to masonry being stronger in compression than in 
tension.) Its elegant shape, more natural and graceful than the straight lines 
and perpendiculars of traditional man-made structures, made it fashionable 
among architects and civil engineers. Its influence can be cited in diverse 
cultures, among which include the Egyptians, Mesopotamians, Romans, (see 
Figure 1.1), Byzantines, French, Chinese, and English. 

A number of investigations on the optimization of arches have been 
conducted over the years. In 1976, Farshad [Ref. 1], using the calculus of 
variations, derived optimality conditions in the form of nonlinear partial dif- 
ferential equations for hinged-hinged arches. An augmented functional, 
comprised of the total potential energy of the system and the objective 
function, appended to the system via Lagrange multipliers, when minimized 
with respect to state variables and with respect to design variables yield the 
system equilibrium equations, and the optimality conditions respectively. 


Three objective functions were imposed: 
— optimal thrust 


— minimum length of arch 


— minimum volume 


The arch span and the loading are specified. The nonlinear system of 
optimality equations were presented but not solved. 

In 1980, Rozvany et al [Ref. 2] considered the problem of arch 
optimization using the Prager-Shield criteria. Here, the arch was in fact a 
funicular frame with beams rigidly interconnected to one another. Only 
statically determinate systems were investigated. The first "arch" with a 
specified span consisted of two inclined beams with a concentrated load along 
the center of symmetry. The second investigation dealt with an "arch" 
consisting of three beam segments, the center segment being horizontal, and 
inclined members from the hinged supports. Two concentrated loads were 
applied at the intersections of the inclined members with the horizontal 
center member. For the single load "arch" it was found that the optimal 
“arch” develops either bending only or axial forces only in the entire structure 
depending on the range of the 4L/d ratio, where L is the span of the structure 
and d is the constant depth of the cross-section. For 4L/d greater than 8, the 
optimal structure has a height half of the span, and there is only axial force 
throughout the structure. For 4L/d less than 8, the optimal structure is a 
straight horizontal beam (i.e., the height is zero), and only bending 
throughout the structure. The width of the beam segments for the optimal 
“arch” varies linearly from the hinged support to the center line. 





Figure 1.1: Pont du Gard of the Nimes Aqueduct 











In a 1980 paper, Lipson et al [Ref. 3] investigated the optimal design of 
arches using the complex method. Both the arch shape and the cross 
sectional dimensions are the design variables for the minimum weight 
structure. Only symmetric arches with constant depth and constant width 
were considered. The section was taken as a thin walled rectangular tube 
with vertical and horizontal wall thicknesses as design variables. The arch 
was approximated with equal length straight beam sections. Thus, each 
beam segment had three design variables, the two wall thicknesses and the 
left end vertical location. In addition, the uniform height and uniform width 
of the rectangular tube were design variables. Side constraints in the form of 
upper and lower bounds were placed on all of the design variables. A 
modified version of the complex method of Box was used as the scheme to 
obtain a "fully-stressed"” optimum design. The shape of the arch was taken as 
a parabola. The optimization algorithm provided the minimum weight 
parabolic arch for a uniform load over a specified span. It was shown that a 
parabolic arch with a rise 0.342 times the span length is the optimal 
parabolic shape for the case of a uniform horizontal load. Deviations within 
10% from this rise have negligible effect on the optimal weight. It was 
further shown that parabolic steel arches will fail due to their own weight at 
span lengths greater than 1,543 ft. For relatively high arches, the maximum 
axial thrust, which occurs at the supports, approaches half of the total 
uniform load. The results of a parametric study of optimal steel arches are 
presented. 

In a 1988 paper, Ang et al [Ref. 4] investigated the optimal shape of an 
arch under bending and axial compression. The cross-section of the arch was 


rectangular with specified constant depth and variable width. With the 


a 








centroidal axis of the arch given by the summation of products of cubic spline 
functions with linear coefficients, the arch axis is permitted to take on any 
smooth shape. The linear coefficients in these products are design variables 
to be determined in the optimization process. The authors considered arches 


under a uniformly distributed horizontal load with three types of boundaries, 
— simply supported-simply supported 
- clamped-clamped 
— clamped-simply supported 


A yield failure constraint was imposed. A new technique for smoothing the 
objective function is presented. It was shown that the optimal shape of the 
arch is a parabola with a rise equal to 0.433 time the span of the arch. No 
other results are presented. It should be noted that the results of [Ref. 3] and 
[Ref. 4], with regard to the ratio of rise to span for an optimal parabolic arch, 
do not agree. 

In the study of arches, it is necessary to determine how they will be 
defined. One prevalent school of thought defines an arch as a curved 
structure, which when supported at both ends and loaded vertically develops 
horizontal reactions. This is apparently intended to eliminate thick walled 
curved beams and straight beams which develop (virtually) no horizontal 
reactions when loaded laterally. 

Others define an arch as a curved beam whose cross-sectional 
dimensions are small relative to its radius of curvature. Hence, the 
centroidal and neutral axes are assumed to coincide. How the structure is 


loaded and supported becomes secondary. This description was chosen to 
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facilitu.te the development of a finite element code capable of generating 
horizontal and vertical displacements and slopes for an arbitrarily loaded 
arch. Without the thin depth assumption, complications arise in the 
calculations of the slopes and displacements because the arch will no longer 
behave as predicted by the beam equilibrium equation: 

(Elv")"= p,(s) (1.1) 
and the bar equilibrium equation: 

(AEu’)’= ~p,(s) (1.2) 
where E represents Young’s Modulus, I the cross-sectional moment of inertia, 
A the cross-sectional area, v the lateral displacement, u the axial 


displacement, p, the axial loading, and Py the lateral loading. 


Once the displacements and slopes are determined, the local stresses 
can be calculated throughout the member using appropriate stress- 
displacement relations. Thus able to determine the stress distribution, the 
arch may be designed to achieve minimum volume (and hence weight) while 
maintaining the developed stresses below some predefined maximum 
allowable stress value. 

The aim of this study is to “optimize” a linear, elastic, isotropic, and 
homogeneous arch under a variety of loadings and end conditions. Although 
these limitations are not physically essential, they were necessary to make 
the investigation tenable given the time constraints of thesis research 
activity. Optimization in this investigation will refer to the variance of the 
cross-sectional geometry to achieve a more uniform stress distribution 
throughout the member. This results in less material used and a more 


efficient structure. VMA Engineering’s Design Optimization Tools [Ref. 5) is 





used to perform the optimization subject to prescribed constraints on the 
design variables as well as stress limitations. The objective function to be 
minimized is the total volume of the arch while maintaining stresses below 
the yield strength of the arch material. Though a rather simplistic model, it 
forms a foundation upon which further research into more complex 


geometries and conditions may be developed. 





II. PROBLEM FORMULATION 


Perhaps the most common optimization in structural mechanics is the 
minimization of an element's weight, subject to a specified loading. Such will 
be the case for this investigation. To make the investigation tenable, the 
problem needed to be narrowed down in its scope. The assumptions and 


approximations made in this study are: 


— The arch is approximated by a series of straight bar/beam elements 
which behave according to the beam equation (1.1) and bar equation 
(1.2). (See Figure 2.1) 


~ The arch material is isotropic, homogeneous, and linearly elastic. 


— The arch’s cross-sectional area will always be of a solid rectangular 
geometry. 


— The arch has a constant circular radius of curvature. 


- The arch “fails” if its internal stresses exceed the yield strength, S,- 


It should be noted that the third and fourth assumptions are not inherent to 
the general optimization problem but rather are imposed to limit the scope of 


this initial investigation. Follow on investigations will relax these 


restrictions. 




















Figure 2.1: Bar Beam Model of the Arch 


Although these suppositions limit the applications for which the optimization 
can be used, they form a foundation upon which further research can be 


based. 
A. THE OBJECTIVE FUNCTION 


The objective of this investigation is the minimization of the structure’s 
weight while maintaining a stress distribution which does not exceed the 
yield strength. Since the weight of the arch is directly proportional to the 
volume of material from which it is made, the objective will be satisfied if the 


total volume of the arch is minimized. That is: 
NEL 
Objective = MIN2 bbl, (2.1) 


where b,, h,, and |, represent the width, height, and length of the i” element, 
respectively, and NEL is the total number of elements. (See Figure 2.1b) 
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Figure 2.2: Arch Dimensions 


With an objective function defined, the next step in the optimization 
process is to impose any necessary (design) constraints upon the system. The 
constraints must be provided to represent the assumptions contained in the 
modeling equations. They should be utilized to avoid undesirable behavior 
such as buckling and yielding. They may also be used to apply any 
limitations on behavior as desired by the designer. The constraints employed 


in this investigation follow. 
B. THE STRENGTH CONSTRAINT 


This study assumes a linearly elastic arch. Therefore, the applied 
loading must not cause the structure to exceed the elastic limit of the 


material from which it is made. Hence, a design constraint which will 
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prevent this mode of "failure by yielding" must be imposed. To fulfill this 
constraint, the internal stresses developed must remain below the yield 
strength of the material. Defining the maximum stress developed in the i** 
element of the arch to be o, and the yield strength of the arch material as S,» 


this constraint becomes 


6; <8, 
or in dimensionless form: 
ofS, -1<0 (2.2) 


C. GEOMETRIC CONSTRAINTS 


To use the beam and bar equations, limits must be imposed upon the 
geometric dimensions of the structure. It is compulsory that the depth and 
width be of at least an order of magnitude smaller than the radius of 
curvature for the beam and bar equations to be applicable. Hence, 

h, < R/10 
or in dimensionless form: 
10h/R, - 1<0 (2.3) 

Similarly, as the width increases, the arch approaches the geometry of 

a shell. To maintain the geometry of an arch, an imposition upon the width 


dimension is also necessary. To avoid shell behavior, a third constraint is 


imposed 
b, < 3h, 
or in dimensionless form: 
b/3h, -1 <0 (2.4) 
ae 











D. SIDE CONSTRAINTS 


Finally, side constraints must be assigned to the dimensions of the 
design variables. For the sake of simplicity, this investigation will only take 
up the variation of the cross-sectional width dimension b;. The side 
constraints must ensure real solutions are obtained, i.e., the arch is a 
physical object and therefore h; and b; cannot be !ess than a realistic finite 
value. The side constraints chosen to reflect these limitations include: 

0.03 in. sh, (2.5) 
0.03 in. <b; (2.6) 

Other constraints could have also been considered. Global buckling 
and local crippling are to name but two. However, the cases to be studied do 
not warrant such a thorough delineation. Therefore, the design and side con- 


straints have been limited to those cited. 
E. OPTIMIZATION SOFTWARE 


With a multitude of preprogrammed optimization routines available, 
the Design Optimization Tools (DOT) software was chosen. Its selection was 
based upon availability, ease of use, and reputation. DOT is a FORTRAN 77 
optimization software package available from VMA Engineering. To perform 
a variety of optimization tasks, DOT uses: 

— The Modified Method of Feasible Directions, 

— Broydon-Fletcher-Goldfarb-Shanno (BFGS) Variable Metric Method, 
— Polynomial Interpolation with bounds, and 

— Sequential Linear Programming (SLP) 
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A user provided “main” program is used to input the variables required by 
DOT. DOT in turn calls a user provided subroutine which defines the 
objective function, design constraints, and design variable side constraints. 
DOT iteratively evaluates the objective function, refining the design variables 
until the optimal solution is obtained. 

The parameters used to calculate the objective function and constraints 
must be known before any optimization can occur. The variables from 
equations (2.1) through (2.6) include: 

— The number of elements used, NEL. 

— The arch radius of curvature, R. 

— The height of the it* element, h.. 

— The length of the i> element, l,. 

— The width of the it> element, b.. 

— The yield strength of the material from which the arch is made, S,: 


~ The stress at the i'* node, o,. 


Of these terms, the number of elements is the choice of the designer. The 
arch radius of curvature and height are constant throughout the span of the 
arch and are defined by the problem. For simplicity, the length of each 
element is made uniform such that: 

1, = @R(NEL) (2.7) 
where 8 represents the subtended arc of the arch. The yield strength is 
determined by the material used to build the arch and the width is the design 


variable to be determined by DOT. 








The stress distribution is not so readily available. However, using the 
beam and bar equations, a finite element scheme can be developed to 
determine the arch’s displacements and slopes due to a given loading. 
Knowing how the displacements and slopes change throughout the arch, the 


stress at a given point can be calculated and the optimization performed. 
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III. STRESS ANALYSIS 


The objective of this optimization is to minimize the total weight 
(volume) of a load bearing arch subject to specified constraints. The strength 
constraint requires that the stress at any point does not exceed the yield 
strength of the arch material. To avoid violating this constraint, the value of 
the stresses at any point must be known. With this requirement, the stress 
development is pursued. 

The normal stress at any point in the arch is composed of normal 
stresses due to bending (bending stress) and normal stresses due to the axial 
forces (axial stress) acting upon the individual elements. Figure (3.1) depicts 
these force interactions. The total normal stress is the algebraic sum of these 
components. 

GO, = 6,+0, (3.1) 
The arch can also develop shear stresses due to shear forces. Due to the 
geometric constraint defined by equation (2.3), the shear stresses turn out to 
be insignificant when compared to the normal stresses. Consequently, the 
shear stresses are ignored.! 

To calculate the normal stress components, we must first determine 
how the elements behave. For a simple straight beam element, the maximum 
normal stress due to bending occurs at the outer fibers and is given by 

6, = Mc/I 


1. See Appendix A for a justification of the omission of the shear stresses. 
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or in terms of the beam equation: 

6, = (Elv")c/I 
which reduces to: 

6, = Ecv" (3.2) 
where c is the distance from the neutral axis to the outer fiber of the beam, 
that is c=h/2. See Figure (3.1b). 





Figure 3.1: Normal Stresses due to Bending and Axial Force 


In the same manner, the normal stresses due to axial behavior can be 
determined. The uniform normal stress due to axial forces, F, acting upon a 
bar can be defined by: 

0, =F/A 
or in terms of the bar equation, 
0, = (AEu’VA 


-—16- 


which reduces to: 
o, = Ew’ (3.3) 
See Figure (3.1c). 
Substituting equations (3.2) and (3.3) into equation (3.1) yields the 
linear equation 
6, = Ecv"+ Eu’ 
or simply 
6, = E(cv" + w) (3.4) 
where v" and wu’ are to be determined from the solution of equations (1.1) and 
(1.2). 

From this development, we see the normal stress is a function of 
Young’s Modulus, the height of the beam, the first derivative of the axial 
displacement, and the second derivative of the lateral displacement. Using 
the Galerkin finite element method, approximate values of u’ and v" can be 
determined. With these values, the stress distribution can be calculated 


using Equation (3.4) and the optimization may then be conducted. 
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IV. FINITE ELEMENT ANALYSIS 


In order to determine the stresses developed for a given loading, the 
values of u’ and v" must be determined. These derivatives can be found by 
solving the beam and bar equations using the Galerkin finite element method 
(FEM). The Galerkin FEM is capable of directly solving systems of linear 
differential equations while preserving their symmetry. 


A. THE BEAM EQUATION DEVELOPMENT 


The beam equation (1.1) is a fourth order linear differential equation 
requiring C! continuity. Therefore, a family of cubic shape functions are 
necessary to maintain function and slope continuity. With this in mind, the 
Galerkin FEM is performed on the beam equation. A finite element method 
is an approximation method which transforms the differential equation of a 
continuous system into a system of linear algebraic equations. The method 
begins by a discretization, that is a partition, of the continuous domain into a 
segmented domain of NEL elements. Thereafter, a three step process takes 
place. 

The first step is to form an approximate solution 9, 

v=v=Qly (4.1) 
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where v is the exact solution in continuous space of the beam equation, v is 
the approximate solution in discrete space, Q™ is the transpose of the column 
vector of cubic shape functions which have the Kronecker delta property, and 
v is the vector of coefficients of lateral displacements and slopes. 
The second step is to form the residual of the approximation where: 
R=£()-p (4.2) 
where £ denotes the system operator, in this case being the beam equation 
such that 
Lv) = [Elv)"]" 
Substituting the beam equation (1.1) and the equation (4.1) into equation 
(4.2) yields: 
R =[E1(QTy)"]’- p,(s) (4.3) 
The third step is to form the Galerkin Equations, 


Jouas = 0 (4.4) 


where Q represents a vector whose values are zero. Substituting equation 


(4.3) into equation (4.4) yields: 


J,ale1Q"y)"T'ds- J,9p,(s) ds = 0 (4.5) 


Performing two successive integrations by parts upon equation (4.5) yields: 


g(E1(Q7v)')', - VEN QTY h+ |,Q"EI(QTv)'ds - J,@p,(s) ds=0 (4.6) 
where B denotes the boundary values of the structure at each end point. 
Since the coefficients of the solution vector are constants, equation (4.6) 


may be rewritten as: 


glE1(Q")"}‘v|, - QEI(QT) yh, + J,Q"EI(QTY'ds y- J,@p,(s) ds=Q (4.7) 
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Using shear given by V=EIv’”’ and moment given by M=ElIv", we define the 
boundary term load vectors V and M as 
Vv =Q[EW Q"'}y |, (4.8a) 
and 
M = QE(Q7)'y |, (4.8b) 


In addition, we define the system stiffness matrix K® as 


B® = |,Q"EI(QTY ds (4.80) 
and the system force vector F° as, 

R= J p9Pp,(s) ds (4.8d) 
and upon substituting equations (4.8) into equation (4.7) we obtain the 
syster: of linear algebraic equations: 

V-M+Ky-F=0 (4.9) 
Moving the applied internal excitation and boundary terms to the right-hand 
side, such that 


K8y=FP+M-V (4.10) 
and defining the load vector of internal and external applied loads as 
FB =Fb4M-V (4.11) 
equation (4.10) simplifies to the linear system: 
Key =F (4.12) 


where the system bending stiffness matrix, K? is constructed from the union 
of the i=1,...,NEL elemental bending stiffness matrices, k>i and the system 
bending force vector F> is formed from the union of the i=1,...,NEL elemental 


bending force vectors f", 
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B. THE BAR EQUATION DEVELOPMENT 


The FEM application to the bar equation is similar to that previously 
conducted with the beam equation. The bar equation, however, is a second 
order linear differential equation requiring C®° continuity. Hence, only a 
family of linear shape functions are necessary to maintain function continuity 
in the FEM development. 

Again, the basic steps of the Galerkin Method are conducted. 

First, the approximate solution 0 is formed: 

u=0=Glu (4.13) 
where u is the exact solution of the bar equation, u is the approximate 
solution, GT is the transposed column of linear shape functions with the 


Kronecker delta property, and u is the vector of coefficients of axial 


displacements. 
Second, the residual is formed: 
R=£(a)+p (4.14) 
where £ pertains to the differential operator of the bar equation, that is, 
L(u) = [AE(u)'T 
Third, the Galerkin equations are formed: 


J ,G(Rods = 0 (4.15) 
Substituting equation (4.13) and the bar equation (1.2) into equation (4.14) 
yields: 
R = [AE(GTu)’] + p, (8) (4.16) 
Substituting equation (4.16) into equation (4.15) yields: 
J,GLAE(GTu)'Tas + J,Gp,(s) ds = 0 (4.17) 
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Performing a single integration by parts upon equation (4.17) yields: 
AEG(G™y) ly - J,G°(AE(GTu) )ds + J,Gp,(s) = 0 (4.18) 
Again, removing the solution vector u of constant coefficients outside of the 
integral yields: 
AE‘G(GT)uh, - JG TAE(GT Jas u + J,Gp,(s) = 0 (4.19) 
Recalling that the axial force F is AEu’, we define the boundary term vector 
Y, 
U =AE'G(G")u |, (4.20A) 
the system stiffness matrix of the bar KA 
KA = J,G°(AE(GT)ds (4.20b) 
and the force vector associated with internal loading p,, as 
Fe = JGp,(s) (4.20c) 
We obtain upon substituting equations (4.20) into equation (4.19), 
U-K{u+F=0 (4.21) 
Taking all the internal excitation F* terms and boundary load terms U to the 
right-hand side of equation (4.21) yields: 
KAy = F*+U (4.22) 
Defining the load vector FA as 
FA= F*+U (4.23) 
and substituting equation (4.23) into equation (4.22) yields the linear system: 
Ku = FA (4.24) 


where 


KA = Uke (4.25) 
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that is, the system axial stiffness matrix, K4 is constructed from the union of 
the elemental axial stiffness matrices, k", and the system axial force vector 


F* is formed from the union of the elemental axial force vectors f'. 


C. THE ELEMENTAL STIFFNESS MATRIX 


In the FEM code, the global (or system) Galerkin FEM equations (4.12) 
and (4.24) are actually constructed from element considerations as follows. 
First, the arch is divided into NEL straight beam-bar elements as illustrated 
in Figure (2.1). The stress analysis program contains a subroutine which 
constructs the elemental beam and bar stiffness matrices, that is, the k>! 
matrices for bending and k* matrices for axial stiffness. The bending and 
axial elemental force vectors, f >! and f*! are also determined within this 
subroutine. Figure (4.1) illustrates the degrees of freedom in which these 


elemental forces (and displacements) act. 
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see node 1 node 1 
1’ lS 2’ 


Bar Element Beam Element 
Figure 4.1: Beam and Bar Element Degrees of Freedom 
The (4x4) ki and (2x2) k*! matrices of the form: 
ky ki ka kid 
key kat kay kza 


bi pbi ,bi bi 
k31 k32 k33 Ky 


bi bi 3 a bi 
ka kl Bla kh 


ai ai 
Ky ky2 


kb = 
= ai ai 
ki k22 
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are combined to form a single (6x6) stiffness matrix, k'. This is accomplished 


by redefining the beam and bar degrees of freedom in the following manner: 


— Redefine the bar local degrees of freedom 1’ and 2’, which refer to the 


axial displacements at each end, as 1’ and 4’ respectively. 














Redefine the beam local degrees of freedom 1’, 2’, 3’, and 4’, which 
refer to the lateral displacement and beam slope at each end as 2’, 3’, 
5’, and 6’ respectively. The redefined degrees of freedom are 


illustrated by Figure (4.2) 


Place the respective components of the beam matrices k>! and bar 


matrices k"' into the elemental stiffness matrix k' where: 


ko «60 kf oC Clo 
o kb Kt og bt xh 
0 KP kB oF KP Kb 
k' o o k# o 0 (4.26) 
0 kP kB oF Kb Kb 


O Ki ke oo kb Kb 


UR 
be 
t 





Beam-bar Element 


Figure 4.2: Beam-Bar Element Degrees of Freedom 
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Defining the elemental displacements and forces in a similar manner, 


the elemental displacement vector becomes: 
(8)? = <8, 8, &, 8, &, a’> (4.27) 
where for the i** element: 
&” = the axial displacement at node 1 
&’ = the lateral displacement at node 1 
&’ = the beam slope at node 1 
5” = the axial displacement at node 2 
& = the lateral displacement at node 2 , 
&Y = the beam slope at node 2 


and are illustrated in Figure (4.2). 
In the same manner, the elemental force vector is redefined as: 
(f')T = <f, of f, i f, i f, Ue f, Me f, v5 (4.28) 


where for the it? element: 
f, ” = the axial force at node 1 


f,” = the lateral force at node 1 
f,” = the moment at node 1 
f, = the axial force at node 2 
f,” = the lateral force at node 2 
f,” = the moment at node 2 
also illustrated in Figure (4.2). 
With these developments, the elemental system of equations for the 


beam-bar element becomes: 
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0 12EI/1} 6EI/13 0 -12EI/1} 6EI/17 le? 
0 6EI/12 4EI/1, 0 -6EI/12 2EI/1, S| jee 
-AE/1, 0 0 AE/I, 0 0 , ‘ " es! 
0 -12EI/1} ~6EI/1? 0 12EI/1} -6EI/1? fz 
C 6EI/17 2EI/1, 0 -6EI/12 4EI/1, . 2 
(4.29) 

or simply 
k's" = £” (4.30) 


Prior to incorporation into the global matrix, a coordinate transformation 


from local to global coordinates is undertaken. 
D. COORDINATE TRANSFORMATION 


Were all the elements of the same orientation with respect to one 
another as it is for a straight beam, a global system of equations could be 
directly constructed. For the arch, however, none of the elements share the 
same orientation. This necessitates the conversion of all elemental 
displacements and forces to a system of global displacements and forces. For 
a reference coordinate system, the horizontal and vertical axes of the arch 
were chosen. (See Figure 4.3) Defining the angle the it element makes with 
the horizontal axis as a,, and the 90° complement of this angle as B., the 
following coordinate relationship between local and global “displacements” 


and "forces" exist, 
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8’ = 8 cos(a,) + 8) cos(B,) 

8 = —8i cos(f,) + 8} cos(a,) 

B= & (4.31) 
8 = 5 cos(a,) + & cos(B;) 

& = Bi cos(f, + & cos(a,) 

& = § 


f, ” =f,‘ cos(a,) + £,' cos(B,) 

f,” = f,! cos(a,) + f,' cos(B,) 

f, ive f, i (4.32) 
f,’ =f, cos(a,) + f;! cos(B,) 

f; vs a cos(@,) + £/ cos(;) 

f,¥ = f,' 






i® element 


Figure 4.3: Displacement & Force Coordinate Transformations 
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Defining [ ‘ as the transformation matrix for the i** element which is 
capable of performing the appropriate coordinate transformations from local 


elemental coordinates to global system coordinates, the matrix becomes: 


cos(a,) cos(B,) 0 0 0 0 
-cos(f,) cos(a@,) 0 0 ) 0 
: r = 0 0 1 0 0 0 
0 0 0 cos(a,) cos(B,) 0 (4.33) 
; 0 0 0 -cos(B,) cos(a,) 0 
0 0 0 0 0 15 
and the relation between local and global “displacements” is 
3; ‘cos(a,) cos(B,) 0 0 ) 0 82 
33 -cos(B,) cos(a) 0 0 0 0) 82 
85 : 0 0 1 0 o Oo gt 
ai 0 0 0 cos(@,) cos(B,) 0 at 
8: 0 0 0 -cos(B,) cos(a,) 0 82 
ai 0 0 oO 0 o 1 82 
(4.34) 


A similar relation between local forces f” and global forces f' exists. The 


notation of equations (4.31) and (4.32) can now be simplified to 


Vv =Lis (4.35) 
fi’ = [ifi (4.36) 
where: 
(S)T = <8,', 8, 5,1, 5,3, 54, 5,'> (4.37) 
(fi)T = <f,', Ee i fe £3 f,i> (4.38) 
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E. THE ELEMENTAL SYSTEM OF EQUATIONS 


Recall from equation (4.30) that k”§” =f". Substituting equations 
(4.35) and (4.36) into equation (4.30) yields 
T'S =Lifi (4.39) 
Multiplying both sides of equation (4.46) by the inverse of the transformation 
matrix, ['', yields 

Ci Big = (Cit igi 
which simplifies to 

(Ciy pp igi = fi (4.40) 
Since [' is an orthogonal matrix, ([ ‘)) = ([ ‘)', and equation (4.41) can be 
rewritten as: 

HBL HS = i (4.41) 
yielding the elemental system of equations transformed to the 
horizontal/vertical coordinate system. Now, the elemental stiffness matrices 
and force vectors are ready for the construction of the global stiffness matrix 
and force vector. The ([ ‘)"k'([ ') term is the elemental stiffness matrix in 


terms of the x- and y-coordinates, and is denoted as k', that is: 


B= (D5)! (C3) (4.42) 
F. THE GLOBAL SYSTEM OF EQUATIONS 


With the elemental system of equations transformed into the global 
(horizontal/vertical) coordinate system, the global system of equations can be 
formed. The system stiffness matrix K is the union of the local transformed 


stiffness matrices for each element, thus 
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where: 





B= (KD) 
and the global (or system) force vector F is obtained by constructing the union 
of the transformed local force vectors f', that is, 

F =uf'i (4.44) 

Then the global system of equations becomes: 
KS =F (4.45) 
Solution of the above system stiffness equations yields the system 
“displacements”. These horizontal, vertical and rotational degrees of freedom 
“displacements” must be transformed back to local axial, lateral and 
rotational ‘digolacements” in order to use the stress equations based upon 
the beam and bar equations. First the global degrees of freedom 6,, 5,, ... , 5,, 
where n=3(NEL+1), are related to the 5; (where i=1,2,...,6) element 
horizontal, vertical, and rotational degrees of freedom for each element. The 


jt degree of freedom for the i** element is given by 
5 = §, (4.46) 

where i=1,2,...,.NEL, j=1,2,3, and k=3(i-1)+j. Then the axial, lateral, and 
rotational "displacements" for the it* element at node 1 are obtained from 
equation (4.31) as 

dy = 5 cos(a,) + & cos(B,) 

3) = —8i cos(B,) + & cos(a,) 

= & 


and likewise for the node 2 end. 








In this manner, the stress at each end of each element can be 
determined. Choosing the greater of the two stresses as the governing stress 
of that element, the optimization analysis can be conducted for the entire 


structure. In this way, the width dimension b, of each element for a 


minimized weight structure is obtained. 
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V. PROGRAM DESCRIPTION AND CAPABILITIES 


Using the previous developments of Chapters II, III, and IV, a 
FORTRAN 77 code was written for execution on the VAX 2000 workstation. 
The program, named ARCH_OPT.FOR, was constructed to be as fully 
interactive with the user as possible to eliminate the need for editing and 
recompiling. The applicability of the code is limited by the assumptions made 
in Chapter II. (i.e., rectangular cross-section, linearly elastic material, etc.) 
As illustrated in Figure (5.1), execution of the program opens and reads an 
input file, ARCH_IN.DAT, which contains information describing the problem 
being investigated. The x-y coordinates of the end points of each element as 
well as the element orientation is determined by the subroutine 
GEOMETRY. The subroutine OPTIMIZATION_TOOL contains the 
parameter OPTDCS, the optimization decision. With OPTDCS=1, DOT is 
called and the weight optimum structure is determined using the provided 
width dimension as the starting point of the optiization process. The stress 
constraint is adhered to based upon the stresses calculated by the FEM 
analysis contained in the subroutine ARCH_STRESS. If no optimization is 
desired, i.e., OPTDCS=0, and the program computes the stress distribution 
based upon the input data, treating the initial geometric parameters as the 
actual design. With the data thus provided, the problem is solved and an 
output file named ARCH_OUT.DAT is created. The output file contains the 


problem parameters, the optimized design variables (width dimensions), and 


the value of the resulting objective function, that is, the minimum volume. 














ARCH_OPT 
GEOMETRY 
:| OPTIMIZATION_TOOL : 
Non-Optimized : . Optimized 


Solution : : Solution 





: : Optimum 
‘ ‘: Design? 


No 


ARCH_OUTPUT 


ARCH_OUT. DAT 


Figure 5.1: ARCH_OPT Program Structure 


To better understand the program’s capabilities, the data fields 
contained in ARCH_IN.DAT need to be discussed. The file is an unformatted 
set of twenty-five numbers separated by commas. This file must be of the 


form: 











ANGLE, RADIUS, YOUNG, YIELD, NEL, METHOD, IPRINT, 
DV1BG, DV1LO, DV1UP, H, CLAN, FX, FY, FM, FA, 
OPTDCS, ITERATE, PRCSN, BX1, BY1, BM1, BX2, BY2, BM2 

Table 5.1 describes each uf these parameters. For further clarification, Figure 


5.2 illustrates how the variables represent the problem and the sign 


conventions used. 
TABLE 5.1: ARCH_IN.DAT FIELD PARAMETERS 


ANGLE A real number from 0 to 180 representing the angle 
subtended by the arch (in degrees). 

RADIUS A real number representing the length of the arch. 
{Dimensions are arbitrary, but they must remain 
consistent for all inputs!} 

YOUNG A real number representing the Young’s Modulus of the 
arch material. 

YIELD A real number representing the yield strength of the 


material used. If a factor of safety is desired, it should 
aby be accounted for and the resultant design strength 
used. 


NEL An integer value from 1 to 32 which denotes the number 


of elements the user wishes to divide the arch for FEM 
evaluation. The prugram is capable of up to 32 elements. 


METHOD An integer from 1 to 2. This is a parameter called by DOT 
to allow the user to select which optimization method is to 
utilized. 

METHOD=1 Modified Method of Feasible Directions 
METi#OD=2 Sequential Linear Programming 


NOTE: If the problem is unconstrained, the BFGS 
algorithm will be used by default [Ref. 5, p. 2-5] 


IPRINT An integer from 0 to 5 used by DOT to control the output 
data from the DOT optimization. See Appendix C for the 
specific outputs 

DV1BG A real number which represents the design variable 1 


(width dimension) best guess. It initializes all element 
width dimension to the best guess value. This establishes 
the optimization starting point. 
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OPTDCS 


ITERATE 





A real number which represents design variable 1’s lower 
limit. (The lower side constraint for the width dimension) 


A real number which represents design variable 1’s upper 
limit. (The upper side constraint for the width dimension) 


A real number which represents the constant height 
(depth) of the arch. 


An integer which represents the number of the node at 
which a concentrated load is to be applied. This number 
must be from 1 to NEL+1. If no concentrated load is 
desired, FX, FY, and FM should be made to equal zero. 


A rea] number which represent the magnitude of a 
concentrated load in the horizontal direction. FX is 
applied at node "CLAN". 


A real number which represent the magnitude of a 
concentrated load in the vertical direction. FY is applied 
at node "CLAN". 


A real number which represent the magnitude of a 
concentrated moment. FM is applied at node "CLAN". 


A real number which represents the magnitude of a 
uniformly distributed lateral load which spans the entire 


length of the arch. 


An integer value which represents the optimization 
"decision’ such that: 


OPTDCS=1 Optimize the dimension of the problem 


OPTDCS=2 Do not optimize the problem. This choice 
will calculate the stress distribution of the 
arch based upon the current problem 
dimensions, assuming the width dimension 
to be constant and equal to DV1BG 


An integer value which represents the number of times 
the resulting "optimized" variables are to be re-entered 
into DOT and the optimization performed again. This 
technique was found to be most useful in refining the 
optimized solutions. 
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PRCSN An integer value from 1 to 2. This parameter allows the 
user to solve the FEM linear system of equations in : 


PRCSNz!1: _ single precision 
PRCSN=2: double precision 

B_.1 An integer value represent the boundary condition of the 
arch’s first nodal point such that if the value is 0, The end 
is free to move in that degree of freedom. If the value is 1, 
the end is fixed in that degree of freedom. 
BX1 = horizontal displacement at point A 
BY1 = vertical displacement at point A 
BM1 = slope of the beam at point A 

B_2 An integer value represent the boundary condition of the 
arch’s last nodal point such that if the value is 0, The end 
is free to move in that degree of freedom. If the value is 1, 
the end is fixed in that degree of freedom. 

BX2 = horizontal displacement at point C 


BY2 = vertical displacement at point C 
BM2 = slope of the beam at point C 


In summary, for a given geometry, loading, and set of boundary 
conditions, the program is able to determine the optimum width dimension of 
each element throughout the length of the arch. This results in an arch of 
minimum volume (weight), capable of supporting the given loading. If 
desired, the optimization process can be bypassed completely. This results in 
the determination of the arch stress distribution based upon the input 
parameters, that is, the stresses associated with the initial design 
dimensions. These factors combine to make ARCH_OPT a useful tool in 


evaluating a variety of arches in engineering applications. 
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Figure 5.2: ARCH_IN.DAT Variable Implementation 





VI. PROGRAM VERIFICATION 


To verify the finite element code used in the optimization investigation, 
several non-optimum problems of beam and arch structures with known 
analytical solutions were solved using the code. The code solution of these 
problems also served the purpose of establishing a relation between the 
number of elements used in the model and the accuracy of the method. With 
a "yardstick" thus provided, "ARCH_OPT.FOR"’s capabilities for accurate 
modeling of beams and arches was assumed. 

The first verification problem was a cantilever beam subjected to a 
concentrated end load, illustrated by Figure (6.1). Gere and Timoshenko 
(Ref. 6, p. 737] give general formulas for the lateral displacements and beam 


slopes for this case, as: 


v = Px®(3L-x)/6EI (6.1a) 
v’ = Px(2]-xV2EI (6.1b) 
Using the parameters: 
P = 1000 lbf h = 3.0 inches 
L = 45 inches E = 30 x 10° psi 
b = 1.5 inches I = bh3/12 


the lateral displacements and beam slopes at the midpoint and free end were 
calculated. ARCH_OPT was run for this beam structure using an angle of 
45.0 x 10°§ radians and a radius of 10° inches to approximate the straight 45 
inch beam length. For a four element approximation, Table (6.1) compares 


the FEM solution to the analytic solution. 
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Figure 8.1: Verification Problem #1 
TABLE 6.1: VERIFICATION PROBLEM #1 SUMMARY OF DISPLACEMENTS 






[ Node 
_savse-oz | S.0008-01 | 
an 5 0. | 2.000800 | —AoTSERo2 3.000E-01 | 
[% Error | | __0.00E+00 | 0.00E+00 | 
ee 5 0. | s.00gee00 | 7.499E-03 9.999E-03 
[% Error | [| 1.33E-02 | _1.00E-02 | 





Max % Error 


where percent error is defined as, 


% Error = 100 x (8,04 — Sep)! Soxact (6.2) 


The stress corresponding to the same points of interest was calculated 
using equation (3.1). Recall 
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which in terms of the beam equation at the i‘ node becomes 


(o,), = (Elv"h/2)/b,h9/12) 


or simply 
(6,); = 6EIv"/b,h? 
From equation (4.1) of the FEM development, v = Q™y, hence 
v" = (QTy)" 
v'(x,) = (QT)"y = (-6/)y; + (4/)y,’ 
Substituting equation (6.4) into equation (6.3) yields 
(o,); = (6EV/b,h?)[-6v,/1? + 4v,1] 


(6.3) 


(6.4) 


(6.5) 


where vy, and v,’ are the lateral dispiacement and slope at x, and are obtained 


from &¥, 8’, & and & of equation (4.31). The stresses of equation (3.1) were 


then compared to those calculated by the ARCH_OPT using equation ( 


6.5). 


The results, summarized in Table (6.2), show a maximum difference of 


5.00x104% between the exact and FEM svlution. 
TABLE 6.2: VERIFICATION PROBLEM #1 SUMMARY OF STRESSES 


Verification Problem #1 


[Node Fixed End_ | Mid Point | __Free End J 
Analytic « | 2.000E+04 | 1.000E+04 |  0.000E+00 | 
FEM o bees 9.999E+03 | 9.712E-06 | 
[% Error | __5.00E-04 | 5.00E-04 | NYA 


. Max % Error aS dee 



















The second verification problem considered was that of a prismatic bar, 


fixed at one end and subjected to an axial concentrated load at the free end as 


shown in Figure (6.2). 
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Figure 6.2: Verification Problem #2 


The parameters used were: 
P = 1000 lbf h = 3.0 inches 
L = 45 inches E = 30x 10° psi 
b = 1.5 inches I = bh3/12 


From the Bar Equation (1.2), we he ve: 
(AEu’y’ = p,(x) 


AEw’ = F(x) (6.6) 
and the normal stress developed in a bar due to axial loading is: 
6, = F(xVA (6.7) 
Substituting equation (6.6) into equation (6.7) yields ; 
6, = AEv/A 
or simply 
6, = Eu’ (6.8) 


Substituting equation (4.18) into equation (6.8) yields the FEM solution, 
6, = E(GTy)’ (6.9) 








From the FEM development, 
(GTuy = (uf) 
which when substituted into equation (6.9) yields the stress at the it* node, 
(o,), = Eu/l (6.10) 
where u,, the axial displacement at x;, is obtained from 5 and 5 of equation 


(4.31). Again using the midpoint and free end as the locations for conducting 
the displacement and stress analysis, the FEM code results are compared to 
those generated by the analytic solutions. These results are listed in Table 
(6.3). They show no difference between the FEM approximation and the 


analytic solution. 
TABLE 6.3: VERIFICATION PROBLEM #2 SUMMARY OF RESULTS 


Verif ication Problem #2 





[Node __|__Fixed End _Mid Point | Free End | 
Analytic 3 0.000E+00 1.667E—04_ 3.333E-04 | 
FEM 3. 0.000E+00 1.667E-04 3.333E-04 | 


[% Error st —sifixed =| (0.00 E+ 00 0.00E+00 | 


[Analytic 5 0.000E+00 ~ 2.222E+02 2.222E+02 | 
FEM 8’ 0.000E+00 2.222E+02 2.222E+02 | 
[% Error) fixed 0.00E+00 | iO. oon | 


“Max % Error 





















The third verification problem chosen was a cantilever beam with a 
concentrated moment at the free end, illustrated by Figure (6.3). Again, from 
Gere and Timoshenko [Ref. 6, p.737], the analytic solution for the 
displacements and slopes of this particular problem are given by the 
equations: 

v = M,x?/2EI (6.11a) 


v= MxEI (6.11b) 


he 








. 
a 
3 
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Figure 6.3: Verification Problem #3 
Using the parameters: 
M, = 10,000 Ibf h = 3.0 inches 
L = 45 inches E = 30 x 10® psi 
b = 1.5 inches I = bh*/12 


the displacements and slopes were determined for the points of interest. 
These results are compared to the 4-element FEM approximations in Table 
(6.4). 

TABLE 6.4: VERIFICATION PROBLEM #3 SUMMARY OF DISPLACEMENTS 


Verification Problem 43 


= Displacement 
ixe 


[Node ss) Fixed End | Mid Point | Free End J 
Analytic 8. 0.000E+00 2.500E-02 1.0000E-01 
FEM 3 0.000E+00 2.500E-02 9.9999E-02 


[1% Error fixed | 0.00E+00 | i100 E-03._| 


Analytic 8 | ae ee 4.444E-03 | 
FEM 8. 0.000E+00 4.444E-03 
[% Error sf fixed 0.00E+00 


E+00 
Max % Error GES03: 

















A comparison of the FEM approximations and analytical solutions 


again show virtually no disparity. These results are presented in Table (6.5). 
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TABLE 6.5: VERIFICATION PROBLEM #3 SUMMARY OF STRESSES 


Verification Problem #3 


[Node Fixed _End Mid Point 
}Analytic o | 4.444E+03 4.444E+03 
| FEM o 4.444E+03 4.444E+03 




















4.444E+03 










Max % Error 





The final verification problem is based upon an example arch problem 


presented by Gere and Timoshenko [Ref. 6, p.616]. They demonstrate how 
the unit-load method can be used to calculate the horizontal displacements of 
the problem illustrated in Figure (6.4). The following formula was obtained: 


8, = PR9/2EI 


y Pp 





7. 


Figure 6.4: Arch Verification Problem #4 
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For the parameters: 
P = 1000 lbf 6 = 90° 
R = 45 inches E = 30 x 10° psi 
b = 1.5 inches I = bh9/12 
h = 3.0 inches 


the horizontal deflection is found to be 0.4500 inches. The horizontal 
deflection from the 4-element FEM approximation is 0.4470 inches, an error 
of 0.66%. 

In all of the verification problems, the percent differences between the 
values obtained by the approximate FEM method and the exact solutions 
were in all cases less than 0.66%. Satisfied that the program was producing 


very good data, the investigation to obtain optimum structures was pursued. 


—~46— 











VII CASE STUDIES 


Recognizing an arch to be but another structural means of transferring 
a load from one point to another, the desire to compare the efficiency of the 
optimized arch to that of a traditional structure drove the first two cases 
studied. Given the problem of transferring the load at B to A, as illustrated 
in Figure (7.1a), numerous structures could be used. For brevity, only the 


frame (7.1b) and arch (7.1c) will be studied. Given the parameters: 


E = 30x10® psi h= 2 inches 
S, = 52,000 psi a = 32 inches 
I = bh?/12 b = 32 inches 


only the width dimension for each case will be allowed to vary. In this way, a 
volume comparison of each structure, hence a measure of the relative 
efficiency of the structure, may be made. 

For the non-optimized frame shown in Figure (7.1b), the base (or width) 
dimension is considered constant throughout. In order to keep the maximum 
stress below the yield strength of the arch material, we have 

(,)max SS, 
or 
Max’! S 52,000 psi (7.1) 
The maximum moment occurs uniformly along the vertical member of the 


frame, hence 


[(2000 Ibf)(32 inX2 in/2)/[b(2 in)3/12] < 52,000 psi 


ah T = 





or simply 
1.846 in. sb (7.2) 
The total volume of the frame is 
Volume = (32 inX2 inX1.846 in) + (32 inX2 inX 1.846 in) 


or 
Volume = 236.3 in® (7.3) 
This volume will be the basis upon which the following optimized 


volumes/weights hence efficiencies will be based. 


2000 libf 2000 lbf 2000 lbf 


[ 
'<___ 


cer eras ensenwcang 


| 


A A 
(a) (b) (ec) 


Figure 7.1: Methods of Transferring a Load 
A. CASE 1: THE OPTIMIZED FRAME 


Given the problem presented in Figure (7.1), an optimized frame of 
equal load bearing capabilities was eought. First, the frame was divided into 
its vertical and horizontal members. The vertical member is subjected to the 


concentrated moment, aP, at C, resulting in a uniform moment and hence 
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uniform stress along the member. Consequently the vertical member has a 
uniform width dimension of 1.846 inches as previously determined. 

The horizontal member is a cantilever beam subjected to a 
concentrated end load. Since the moment along BC varies linearly from 0 at 


B to aP at C, for a constant Srmax= Sy? the width must also vary linearly. Thus 


the width varies linearly from 0 inches at the right end to 1.846 inches at C. 
The volume for this frame is: 
Volume = (32 inX2 in)(1.848 in) + (.5X32 inX2 inX1.848 in) 
which is: 
Volume = 177.4 in® (7.4) 

The volume of the horizontal member, BC, was then optimized using a 
4-element, 8-element, and 12-element discretization. Table (7.1) illustrates 
the optimized volumes of each of these solutions. The percent difference 
between the 4-element and 8-element solution was found to be less than 
11.5% and that for the 8-element and 12-element solution to be less than 3%. 
In the interest of solving many cases, it was decided to solve all future case 
studies with a 12-element discretization, treating the 12-element model as 
producing grid independent results. 

The optimized results for the horizontal member, given in Table (7.1) 
show a total member volume of 64.63 in. The volume of the vertical member 
remains the same as for the non-optimized structure, hence the total volume 
of the optimized frame is 

Volume = (32 inX2 inX 1.848 in) + 64.63 in? 
= 182.9 in® (7.5) 
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This represents 22.6% less volume, hence 22.6% less weight than the non- 


optimized frame. 
TABLE 7.1: CASE 1 SUMMARY OF RESULTS 


honed FE PP PP id Ba 
inches inches inches cubic in. pai 


















































































1 2.000E+00 2.667E+00 1 ‘B4BE100 9.856E+00 §.19E+04 
2 2.000E+00 2.667E+00 1.6965E+00 9.040E+00 ; 5.19E+04 
3 2.000E+00 2.667E+00 1.543E+00 8.229E+00 3 5.1BE+04 
4 2.000E+00 2.687£+400 1.389E+00 | 7.40B8E+00 4 5§.1BE+04 
5 2.000E+00 2.667E+00 1.236E+00 6.592E+00 5 §.1BE+04 
6 2.000E+00 2.667E+00 1.082E+00 5.771E+00 8 5.18E+04 
7 2.000E+00 2.687E+00 9.313E-01 4.987E+00 7 5.17E+04 
B 2.000E+00 2.667E+00 7.7B1E-O1 4.134E+00 8 §.15E404 
9 2.000E+00 2.667E+00 6.316E-01 3.369E+00 9 §.16E+04 

2.000E+00 2.667E+00 4676E-01 2.494E+00 6.07E+04 

2.000E+00 2.667E+00 3.196E-01 1.705E+00 5.138404 





5.07E+04 
4.00E+04 


2.667E+00 : 
12-element Z Volume. 228: ee 
B-element Z Volume: 6 '849E40 1 
4-element Zz Volume: 7.385E+01 


2.000E+00 


B. CASE 2: THE OPTIMIZED CANTILEVER ARCH 


Given the circumstances and parameters of Figure (7.1), a cantilever 
circular arch (Fig. 7.1c) was employed to perform the same function, 
transferring the given load at point B to point A. The resulting dimensions 
and stresses of the optimization are presented in Table (7.2). These results 
illustrate what one would have expected, the width dimension is reduced 
until the local stress approaches the yield strength of the material. The total 
volume of the arch, 128.3 in® is 46.3% less than the non-optimized frame and 
29.9% less than the optimized frame. In moving structures where higher 


weights mean higier operating costs, savings such as these can become 


significant. 
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TABLE 7.2: CASE 2 SUMMARY OF RESULTS 


oA eo el 
inches nches inch Bat na sary 























































1 4.186E£+00 1.5618+01 5.20E+04 
2 4.1B6E+00 1.542E+01 ; 5.20E+04 
3 4.1B6E+00 1.503E+01 3 5.20E+04 
4 4.186E+00 1.438E+01 4 5.20E+04 
5 4.186E+00 1.350E+01 5 5.19E+04 
6 4.1B6E+00 1.237E+01 6 5.19E+04 
7 4.1B6E+00 1.109E+01 7 5.16E+04 
8 4.186E+00 9.762E+00 8 6.16E+04 
] 4.186E+00 7.799E+00 9 5.05E+04 

4.186E+00 6.966E+00 5.19E+04 










4.85E+04 
3.47E+04 
1.79E+02 


5.168E-01 | 4.327640 
3.856E-01 | 3.061E+00 | 
EZ Volume: 2:t:288R362: 


2.000E+00 
2.000E+00 


4.186E+00 
4.186E+00 






C. CASE 3: THE OPTIMIZED CANTILEVER ARCH WITH AXIAL 
END LOAD 


To appreciate how the stress distribution follows the moment 
distribution of the member, the problem of Case 2 was modified keeping the 
same parameters as outlined previously, by changing the direction of the load 
as shown in Figure (7.2). 


me em 








2000 lbf 





Figure 7.2: Case 3 Problem Geometry 


For this case, the moment at any point is given by the expression 

M = PR(1-sin(®)) (7.6) 

as opposed to the moment distribution of Case 2 where the moment along the 
length of the arch is 

M = PRcos(@) (7.7) 

Assuming a unit moment, that is PR=1, then from Figure (7.3), one may see 

that the area under Case 3’s moment curve is significantly less than the area 

under Case 2’s moment curve. One would expect this to correspond with the 


need for less material due to less applied force (in this case, moment). 
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Figure 7.3: Moment Diagrams of Case 2 and Case 3 


The resulting optimized arch of this case has the dimensions outlined 
in Table (7.3) and a noticeably smaller volume of 81.39 in? than the 128.3 in? 
of Case 2. In fact, the reduction in volume is the same as the reduction in 
areas of the moment diagrams. This result collaborates with the previous 
observation. Hence one may see that the normal stress is predominantly due 
to bending and essentially follows the moment distribution of the structure in 


question. 
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TABLE 7.8: CASE 3 SUMMARY OF RESULTS 


mie | ieee | ame Tes PL 
nche neh inches! suit ie ips 














2.000E+00 
2.000E+00 
2.000E+00 
2.000E+00 











4.186E+00 
4.186E+00 
4.186E+00 
4.1B6E+00 








2.611E-01 
2.000E-01 
2.000E-01 



















2.186E+00 
1.674E+00 
1.674E+00 


1 2. OOF +00 4.186E+00 1.646E+00 1.545E+01 
2 2.000E:00 | 4.1686E+00 1.687E+00 1.412E+01 
3 2.000E+00 | 4.186E+00 1.435E+00 1.201E+01 
4 2.000E+00 | 4.186E+00 | 1.196E+00 1.001E+01 
5 2.000E+00 | 4.186E+00 | 9.703E-01 | 8.123E+00 
6 2.000E+00 | 4.186E+00 | 7.574E-01 | 6.341E+00 
? 2.000E+00 | 4.186E+00 | 5.676E-01 | 4.752E+00 
8 2.000E+00 | 4.166E+00 | 4.009E-01 | 3.356E+00 
9 
10 





D. CASE 4: THE CANTILEVER ARCH UNDER A DISTRIBUTED 
LOAD 


This case involves applying a uniformly distributed load acting radially 
inward on a cantilever arched segment as pictured in Figure (7.4) 


- 


p(s) = 100 lbf/in 





Figure 7.4: Case 4 Problem Geometry 








where: 


E = 330x108 psi h = 2 inches 
S, = 52,000 psi R = 32 inches 
I = bh?/12 © = 90 degrees 


The results of the optimization contained in Table (7.4) illustrate how 
the critical constraints shift from the maximum allowable stress to the 
minimum allowable width dimension. Even though the moment in the 
structure has diminished and the stress no longer approaches the yield 
strength of the material, the geometric constraint prevents the cross-section 
from becoming so thin that the bar/beam assumption is no longer valid. The 
variation in the width dimensions appears to be almost logarithmic alluding 
to the complexities involved with arched segments under uniformly 
distributed loads. Had the arch not been optimized and assuming the 
maximum width dimension of Table (7.4) to represent the uniform width 
dimension of the non-optimum arch, such that 

Volume = [(R,)? — (R,)?(@ radiansXb,,, ,/2 (7.8) 
then the total volume of the non-optimum arch would have been 
Volume = 209.1 in.3 
Hence, the optimized volume of 101.6 in* is 51% smaller than that of the 


non-optimized arch. 


S165 














TABLE 7.4: CASE 4 SUMMARY OF RESULTS 


oe a 
nches inches inches ubic in. 


abr 


































£ 950E+00 
3.266E+00 
1.857E+00 
1.674E+00 
= S7AE 00. 


6.812E-01 
3.901E-01 
2.21BE-O0t 
2.000E-01 


2.000E+00 
2.000E+00 
2.000E+00 
2.000E+00 
2.000E+00 







1 2.000E+00 2.080E+00 1.741E+01 
2 2.000E+00 1.930E+00 1.616E+01 
3 2.000E+00 1.762E+00 1.475E+01 
4 2.000E+00 1.552E+00 1.299E+01 
5 2.000E+00 1.315E+00 1.101E+01 
6 2.000E+00 1.067E+00 | 8.933E+00 
7 2.000E+00 B8.213E-01 , 6.876E+00 
8 

9 






4.186E+00 
4.186E+00 





E. CASE 5: THE SIMPLY SUPPORTED ARCH 


This case involved the optimization of a simply supported arch 
subjected to a lateral load at the midpoint, illustrated in Figure (7.5). The 
arch is a first order statically indeterminate structure subject to the following 


parameters: 
E = 30x10 psi h= 2inches 
S, = 52,000 psi R = 32 inches 
I = bh9/12 © = 180 degrees 
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Figure 7.5: Case 5 Problem Geomctry 


In order to take advantage of symmetry, the problem was divided along the 


axis of symmetry with the following boundary conditions imposed on the 


symmetry end, 
(EIv")’ = P/2 
ElIv’ =0 
u=0 


The results of the optimization are summarized in Table (7.5). Here 
again, the weight savings of the optimized arch of 263.2 in? over that of the 
non-optimized arch, as defined by equation (7.8), of 628.1 in? is approximately 
58%. 


7 ae 





TABLE 7.5: CASE 5 SUMMARY OF RESULTS 


oe a ial 
inches inches ncheg pees a abr 









































































1 2.000E+00 {| 4.186E+00 | 5. 1S0E-01 4.295E+00 
2 2.000E+00 | 4.186E+00 | 9.069E-01 | 7.592E+00 
3 2.000E+00 | 4.186E+00 1.117E+00 | 9.851E+00 
4 2.000E+00 | 4.186E+00 1.478E+00 1.23BE+01 
5 2.000E+00 | 4.186E+00 1.196E+00 1.001E+01 
6 2.000E+00 | 4.186E-00 1.127E+00 | 9.435E+00 
7 2.000E+00 | 4.186E+00 1.531E+00 | 1.282E+01 
8 2.000E+00 | 4.186E£+00 | 5.765E-01 | 4826E+00 
9 2.000E+00 | 4.186E+00 | 5.935E-01 | 4.969E+00 

2.000E+00 | 4.1B6E+00 1.346E+00 1.127E+01 





2.000E+00 
2.000E+00 


4.186E+00 
4.186E+00 


2.206E+00 | 1.847E+01 


: &. 616E+01 








F. CASE 6: THE FIXED-FIXED ARCH 


To determine the effect of additional redundancy, the next case 
involved the optimization of a fixed-fixed arch subjected to a lateral load at 
the midpoint, illustrated in Figure (7.6) and subject to the following 


parameters: 
E = 30x10® psi h= 2 inches 
S, = 52,000 psi R = 32 inches 
I = bh*/12 © = 180 degrees 
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Figure 7.6: Case 6 Problem Geometry 


The fixed supports add two additional redundant moments at the supports 
when compared to the previous simply-supported, simply-supported arch of 
case §. Again taking advantage of symmetry, the problem was divided along 


the axis of symmetry with the following boundary conditions imposed on the 


symmetry end, 
(EIv") = P/2 
Elv’ =0 
u=0 


The resulting optimized arch has a total volume of 210.3 in? as 
summarized in Table (7.6). The non-optimum arch has a volume of 533.0 in? 
assuming a constant arch width corresponding to the optimized arch’s 


maximum width of 2.651 inches. Here, the weight savings of the optimized 
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arch over that of the non-optimized arch is approximately 60%. It is also 
worth particular note that though the loading and geometry of Cases 5 and 6 
are the same, differing only with regards to the boundary conditions imposed, 
Case 6 is 20% lighter. This is because a fixed-fixed structure is more 
statically indeterminate than the simply supported structure. Hence we see 
here that the more redundant structure is the more efficient member. This 
illustrates one of the reasons why fixed-fixed structures are preferred in 
construction. 


TABLE 7.6: CASE 6 SUMMARY OF RESULTS 


aad a a 
inches! inches inches cubic in. pal 















































Yr Volume: 
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2.219E4+01 















1 2.000E+00 4.186E+00 1.727E+00 1.44B6E+01 5.20E+04 
2 2.000E+00 4.1868E+00 | 9.458E-01 7.918E+00 : 5.00E+04 
3 2.000E+00 4.186E+00 2.719E-O1 2.276E+00 3 4.62E+04 
4 2.000E+00 4.186E4+00 | 6.504E-01 5.445E+00 4 5.20E+04 
5 2.000E+00 4.186E+00 1.091E+00 | 9.133E+00 5 4.946404 
6 2.000E+00 4.186E+00 | 7.899E-01 8.613E+00 8 5.20E+04 
7 2.000E+00 4.186E+00 1.062E+00 8.891E+00 7 5.18E+04 
6 2.000E+00 4.1B86E+00 | 7.539E-01 6.311E+00 6 4.21E+04 
9 2.000E+00 4.186E+00 2.585E-01 2.164E+00 9 5.20E+04 

2.000E+00 4.186E+00 | 9.328E-01 7. BOBE+00 10 } 6.138404 

2.000E+00 4.186E+00 1.753E+00 1.46BE+01 5.24E+04 

2.000E+00 4.186E+00 | 26616+00 5.19E+04 














VIII CONCLUSIONS 


The conclusions of this study are: 


the 


The stress analysis based upon the bar/beam model yielded good 
results with percent deviations from known analytic solutions 
ranging between 0.1 and 1.5%. Hence, the bar/beam element model 


is a viable technique in the approximation of arch structures. 


The DOT optimization software was able to utilize the bar/beam 
modeled stress analysis to efficiently determine weight optimum 


arch structures. 


The optimization demonstrates how structures which are more 
statically indeterminate (redundant) are likewise more efficient than 


identical structures under identical loading. 


The weight optimization of a structure is available and effective for 


all types of problem boundary conditions. 


It should be noted once again that this is an initial investigation into 


weight optimization of arches. Numerous opportunities exist for the 


expansion of the basic assumptions made in this study. This investigation 


only considered the optimization of structures of the form: 


I=kA 
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The general form of this type of optimization is such that 

I=kA® 
where n = 1, 2, or 3 depending upon which cross-sectional dimension(s) is 
defined as the design variable. Some of the possibilities for future research 


include: 


- Allowing only the cross-sectional height dimension, h, vary while 





holding all other parameters constant. (n=3) 


- Allowing both the height and width dimensions to vary 
proportionally, while holding all other parameters constant. (n=2) 


— Allowing the radius of curvature of the arch (its center-line shape) to 
vary. 
— Optimizing the arch using engineering cross-sections such as box 


beams, I-beams, circular cross-section, etc. 


— Incorporating additional constraints (such as arch maximum height 
limitations, buckling constraints, crippling constraints, etc,) in order 


to expand the model; thereby enabling the model to solve a greater 


variety of problems. 








APPENDIX A 


JUSTIFICATION FOR OMITTING SHEAR STRESSES 


The shear stress distribution through a beam of rectangular cross- 
section has a parabolic distribution along the height of the member. The 
maximum shear stress, located at the neutral axis of the beam, is 

Tmax = 1.5V/A (A.1) 
where t,,,, is the maximum shear stress, V is the shear force, and A is the 
cross-sectional area of the beam. [Ref. 6, p. 229] 

The normal stress due to bending is given by the equation 

6, = Mc/l (A.2) 
where o,, is the maximum normal stress, M is the bending moment, and I is 
the cross-sectional moment of inertia which for this case is bh?/12 where b 
and h are the width and height respectively of the cross-section. 

Redefining the normal stress in terms of the cross-sectional dimensions 
yields 

6, = M(h/2)(bh9/12) 
or 
0, = 6M/hA (A.3) 

The ratio of the maximum shear stress to tue normal stress due to 

bending, is denoted by r and given by the expression: 
r=1.,,/5, (A.4) 








Substituting equations (A.1) and (A.3) into equation (A.4) yields 
r = (1.5V/AV(6M/hA) 
or 
r= Vh/4M (A.5) 
For the cases investigated in this study, the maximum value r can attain is 
when the loading is that of a uniformly distributed load, Py. Then, where: 
V=p,L (A.6) 
M= p,L*/2 (A.7) 
which upon substitution into equation (A.8) yields 
r= (pyL)h/4(p,L?/2) 
which simplifies to 
r=h/2L (A.8) 
The use of the beam equation requires the length of the beam to be at a 
minimum ten times the height, that is: 
L2 10h (A.9) 
To maximize the value of r, let L equal 10h, the minimum allowable length. 
Cabstituting this value of L into equation (A.8) yields 
r $ h/2(10h) 
or simply 
rs 1/20 (A.10) 
Hence, the maximum shear stress accounts for less than five percent of the 
bending stress developed in the structure. Five percent is high considering 
this analysis over-assumed the value of the shear stress by assigning the 
maximum shear stress to the entire cross-section of the beam. Moreover, at 


the outermost fibers where o, is a Maximum, the shear stress is zero. 











Therefore, under the circumstances of this study, the addition of shear 


stresses was deemed to be unwarranted. 
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APPENDIX B 


DOT USERS MANUEL, SECTION 2.1 


DOT WITH APPLICATION PROGRAMS 





2.1 CALLING STATEMENT 


DOT is invoked by the following FORTRAN calling statement 
in the user’s program: 


CALL DOT ( INFO, METHOD, IPRINT, NDV, NCON, X, 
* XL, XU, OBJ, MINMAX, G, RPRM, IPRM, WK, NRWK, 
* IWK, NRIWK) 











APPENDIX C 


DOT USERS MANUEL, SECTION 2.2 





2.2 PARAMETERS IN THE CALLING 
STATEMENT 


Table 2-1 lists the parameters in the calling statement to DOT. 
Where arrays are defined, the required dimension size is given 
as the array argument. These are minimum dimensions. The 
arrays can be dimensioned larger than this to allow for program 
expansion. 


TABLE 2-1: F.\ RAMETERS IN THE DOT ARGUMENT LIST 


PARAMETER 


Information parameter. Before calling DOT 
the first time, set INFO=0. When control 
returns from DOT to the calling program, 
INFO will normally have a value of 0 or 1. 

tf INFO= 0, the optimization is complete (or 
terminated with an error message). If 
INFO= 1, the user must evaluate the objec- 
tive, OBJ, and constraint functions, G(1), 
iz 1,NCON, and call DOT again. A third 
possibility, INFO= 2, exists also. In this 
case, the user must provide gradient in- 
formation. This is an advanced feature 
and is described in Section 3.2. 


Optimization method to be used. 

METHOD = 0 or 1 means use the modified 
method of feasible directions. 

METHOD = 2 means use the sequential 
linear programming method. If the problem 
in unconstrained (NCON = 0) the BFGS al- 
gorithm will be used, regardless of the 
value of the parameter METHOD. 
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DOT WITH APPLICATION PROGRAMS 


IPRINT Print control parameter. 

IPRINT = 0 no output. 

IPRINT = 1 internal parameters, initia! 
information and results. 

IPRINT = 2 same plus objective function 
and X-vector at each feration. 

IPRINT = 3 same plus G-vector and 
critical constraint numbers. 

IPRINT = 4 same plus gradients. 

IPRINT = 5 same plus search direction. 

NOTE: The IPRM Array contains additional 

print options. 


Number of decision/design variables con- 
tained in vector X. NDV is the same as N 
in the mathematical problem statement 
given in Section 1.7. 


Number of constraint values contained in 
array G. NCON is the same as M in the 
mathematical problem statement given in 
Section 1.7. NCON=0 is allowed. 


| Vector containing the design variables. On 
the first call to DOT, this is the user’s best 
guess of the design. On the final return 
from DOT (INFOx0 is retumed), the vector 
X contains the optimum design. 


XL(NDV) Array containing lower bounds on the 
design variables, X. If no lower bounds are 
imposed on one or more of the design vari- 
ables, the corresponding component(s) of 
XL must be Set to a large negative number, 
say -1.0E+15. Be sure it's -1.0E+15 and 
not -1.0E-15 (+15, not -15). 
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Array containing upper bounds on the 
design variables, X. ff no upper bounds 

are imposed on one or more of the design 
variables, the corresponding component(s) 
of XU must be Set to a large positive num- 
ber, say 1.0 E+15. 


XU(NDV) 











Value of the objective function correspond- 
ing to the current values of the design vari- 
ables contained in X. On the first call to 
DOT, OBJ need not be defined. DOT will 
return a value of INFO=1 to indicate that 
the user must evaluate OBJ and cali DOT 
again. Subsequently, any time a value of 
INFO=1 is retumed from DOT, the objec- 
tive, OBJ, must be evaluated for the cur- 
rent design and DOT must be called again. 
OBJ has the same meaning as F(X) in the 
mathematical problem statement given in 
Section 1.7. 















Integer parameter specifying whether the 
minimum (MINMAX«0,-1) or maximum 

(MINMAX=1) of the objective junction is to 
be found. 








G(NCON) Array containing the NCON inequality con- 
straint values corresponding to the current 
design contained in X. On the first call to 
DOT, the constraint values need not be 
defined. On retum from DOT, if INFO=1, 
the constraints must be evaluated for the 
current X and DOT must be called again. 

i ff NCON=0, array G must be dimensioned 
to 1 or larger, but no constraint values 
need to be provided. 
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IWK(NRIWK) 


2-8 2.2 PARAMETERS IN THE CALLING STATEMENT 





Array containing the real (floating point 
numbers) control parameters. Initialize the 
entire array to 0.0 to use all defautt values. 
If you use other values than the defaults, 
Set the corresponding entries to the 
desired values. Section 3.1 describes how 
to change the value of these parameters. 













Array containing the integer contro! 
parameters. As with the RPRM array, set 
the array to zero to use the default values, 
or set the proper entries to the desired 
values. Section 3.1 describes how to 

change the value of these parameters. 










User provided work array for real (floating 
point) variables. Array WK is used to store 
imernal scalar variables and arrays used 
by DOT. If the user has not provided 
enough storage, DOT will print the ap- 
propriate message and terminate the op- 
timization. 










Dimensioned size of work array WK. 
NRWK should be set quite large, starting at 
about 500 for a small problem. tf NRWK 
has been given too small a value, an error 
message will be printed and the optimiza- 
tion will be terminated. 










User provided work array for integer (fixed 
point) variables. Array IWK is used to store 
imternal scalar variables and arrays used 
by DOT. If the user has not provided 
enough storage, DOT will print the ap- 
propriate message and terminate the op- 
timization. 
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Dimensioned size of work array |WK. A 
good estimate is 300 for a small problem. 
Increase the size of NRIWK as the problem 


grows larger. If NRIWK is too small, an 
error message will be printed and the op- 
timization will be terminated. 





Note: The minimum required values of NRWK and NRIWK are 
defined as follows (The dimensions may be larger than this): 


N1 = NCON +NDV 

N2 = 2*NDV 

N3 = 10*NDV 

N4 = MIN (N1,N2) 

N5S=1 

IF NCON = 0, N5=0 

NCOLA = MAX (N3,N4) 
NGMAX = MIN (NCON,NCOLA) 
NRB = MIN (N1,NCOLA+1) 

IF NCON = 0, NRB= 1 


NRWK = NDV*(10+NCOLA) + 5*NCON + NCOLA+ NRB**2 + 
MAX(NDV,NCOLA)+2*NS5*NRB + 40 


IF METHOD > 1, NRWK = NRWK + NGMAX + 
NDV*(3+NCOLA) 


NRIWK = NDV + NCON + NGMAX + 71 
IF METHOD > 1, NRIWK = NRIWK + NGMAX 


2.2 PARAMETERS IN THE CALLING STATEMENT 2-9 
Prey fs BS 
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A program called DTSTOR is provided with DOT. If you com- 
pile and run this program interactively, the minimum required 
values of NRWK and NRIWK ure calculated for you. See Ap- 
pendix B for more information «<n this option. 
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APPENDIX D 


VERIFICATION AND CASE STUDY OUTPUT 


VERIFICATION #1 


OPTIMIZATION SOLUTION 


A) Problem Parameters: 


Arch Angle : 0.003 Youngs Modulus: 30000000.0 
Arch Radius: 1000000.000 Yield Strength: 52000.0 
Arch Height: 3.000 No of Elements: 4 
B) Derived Constants: 
No of System Nodal Points... 5 
No of Degrees of Freedom.... 15 
Length per Element.......... 11.2500 
Number of Iterations........ 0 
Cc) Structure Loading: 
FR os sic A 3: 6ie wiaheire Sderaid oe eee weetens 1000.0000 
BEY facia fevatla Taal aie idre casing ey ecg tele 0.0000 
EM geo) ats tor bya a ac bre co's tls tad see Sus evans 0.0000 
BA i625 soos oc 5h Sire nose ee a emacs eb wee 0.0000 
D) Elemental Dimensions and Stress Distribution: 
Element Height Base Length Volume 
1 3.00000 1.50000 11.24996 50.62481 
2 3.00000 1.50000 11.24996 50.62481 
3 3.00000 1.50000 11.24996 50.62481 
4 3.00000 1.50000 11.24996 50.62481 
E) Objective Function: 
Total structure Volume: 202.499222 
Node Stress 
1 19999.90 
2 14999.94 
3 9999.967 
4 4999.973 
5 9.7119728E-06 
F) Boundary Conditions: 
Node X-Displ Y-Displ Slope 
1 1 1 1 
5 0 0 0 
G) Solution Vector: 
Node X-Displ Y-Displ Slope 
1 0.000000E+00 0.000000E+00 0.000000E+00 
2 0.257809E-01 0.112328E-08 -0.437496E-02 
3 0.937488E-01 0.409060E-08 -0.749993E-02 
4 0.189841E+00 0.828730E-08 -0.937491E-02 
5S 0.299996E+00 0.130987E-07 -0.999991E-02 


af: 











VERIFICATION #2 


OPTIMIZATION SOLUTION 


A) Problem Parameters: 


Arch Angle : 0.003 Youngs Modulus: 30000000.0 

Arch Radius: 1000000.000 Yield Strength: 52000.0 

Arch Height: 3.000 No of Elements: 4 
B) Derived Constants: 

No of System Nodal Points... 5 

No of Degrees of Freedom.... 15 

Length per Element........-.- 11.2500 

Number of Iterations........ 0 





C) Structure Loading: 


Bia er le Boe ea ala Sense Ghee as 0.0000 
EY ccd sa: greens ewe le Rneete Be. BEES 1000.0000 . 
BM cb deed sie Suace sh ese Bee ed Shas 0.0000 
BA bo Sica ek Si ib Sela e oie wre Saas S 0.0000 
D) Elemental Dimensions and Stress Distribution: 
Element Height Base Length Volume 
1 3.00000 1.50000 11.24996 50.62481 
2 3.00000 1.50000 11.24996 50.62481 
3 3.00000 1.50000 11.24996 50.62481 
4 3.00000 1.50000 11.24996 50.62481 


E) Objective Function: 
Total structure Volume: 202.499222 


Node Stress 
1 222.2231 
2 222.2229 
3 222.2227 
4 222.2224 
bs) 222.2222 

F) Boundary Conditions: 

Node X-Displ Y-Displ Slope 
1 1 1 1 
5 0 0 0 


G) Solution Vector: 
Node X-Displ Y-Displ Slope 


1 0.000000E+00 0.000000E+00 0.000000E+00 
2 0.112328E-08 0.833330E-04 -0.191236E-09 
3 0.409060E-08 0.166666E-03 -0.327832E-09 . 
4 0.828730E-08 0.249999E-03 -0.409790E-09 
5 0.130987E-07 0.333332£-03 -0.437110E-09 








VERIFICATION #3 


OPTIMIZATION SOLUTION 


A) Problem Parameters: 


Arch Angle : 0.003 Youngs Modulus: 30000000.0 
Arch Radius: 1000000.000 Yield Strength: 52000.0 
Arch Height: 2.000 No of Elements: 4 
B) Derived Constants: 
No of System Nodal Points... 5 
No of Degrees of Freedom.... 15 
Length per Element.......... 11.2500 
Number of Iterations........ 0 
C) Structure Loading: 
EX 6.05, ise) bits Ses ars Loe! ois Vote eerie Suatele where 0.0000 
BEY, 0 15e: Seeuaiar sets cbt ay ose ao rahe oud veloss 0.0000 
BM oo sane ane eines "oy ites oe, So eek ee eas 10000.0000 
BA 56 od stave te orate ecoleresatere essence 0.0000 
D) Elemental Dimensions and Stress Distribution: 
Element Height Base Length Volume 
1 3.00000 1.50000 11.24996 50.62481 
2 3.00000 1.50000 11.24996 50.62481 
3 3.00000 1.50000 11.24996 50.62481 
4 3.00000 1.50000 11.24996 50.62481 


E) Objective Function: 
Total structure Volume: 202.499222 


Node Stress 
1 4444.438 
2 4444.438 
3 4444.441 
4 4444.441 , 
5 4444.441 

F) Boundary Conditions: 

Node X-Displ Y-Displ Slope 
1 1 1 1 
5 0 0 0 


G) Solution Vector: 
Node X-Displ Y-Displ Slope 
1 0.000000E+00 0.000000E+00 0.000000E+00 
2 -0.624994E-02 -0.273194E-09 0.111110E-02 
3 -0.249998E-01 -0.109277E-08 0.222221E-02 
4 -0.562495E-01 -0.245874E-08 0.333332E-02 
S -0.999991E-01 -0.437110E-08 0.444442E-02 
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VERIFICATION #4 


OPTIMIZATION SOLUTION 


Arch Angle : 90.000 Youngs Modulus: 30000000.0 
Arch Radius: 45.000 Yield Strength: 52000.0 
Arch Height: 3.000 No of Elements: 4 
B) Derived Constants: 
No of System Nodal Points... 5 
No of Degrees of Freedom.... 15 
Length per Element.......... 17.5581 
Number of Iterations........ 
C) Structure Loading: 
EX easton acelerg aera tens eniie veieielie  tiectveve 0.0000 
BY tiers Sek eie beets eae oe Sa ees 1000.0000 
BM iscik SGihiale Sa ietele So Se Se we eee SS 0.0000 
BA en 838 ob pen 6 io cate or eys 6 Gig eet owes 0.0000 
D) Elemental Dimensions and Stress Distribution: 
Element Height Base Length Volume 
1 3.00000 1.50000 17.55813 79.01159 
2 3.00000 1.50000 17.55813 79.01159 
3 3.00000 1.50000 17.55813 79.01159 
4 3.00000 1.50000 17.55813 79.01159 
E) Objective Function: 
Total structure Volume: 316.046356 
Node Stress 
1 20217.61 
2 18601.18 
3 14266.01 
4 7777 .182 
5 43.39704 
F) Boundary Conditions: 
Node X-Displ Y-Displ Slope 
1 1 1 1 
5 0 0 0 
G) Solution Vector: 
Node X-Displ Y-Displ Slope 
1 0.00C900E+00 0.000000E+00 0.000000E+00 
2 -0.654617E-01 0.131512E-01 0.750656E-02 
3 -0.223502E+00 0.118880E+00 0.138705E-01 
4 -0.381543E+00 0.355536E+00 0.181227E-01 
5 -0.447005E+00 0.684770E+00 0.196159E-01 





OPTIMIZATION #1 


OPTIMIZATION SOLUTION 


A) Problem Parameters: 


Arch Angle : 0.002 Youngs Modulus: 
Arch Radius: 1000000.000 Yield Strength: 
Arch Height: 2.000 No of Elements: 
B) Derived Constants: 
No of System Nodal Points... 5 
No of Degrees of Freedom.... 15 
Length per Element.......... 8.0000 
Number of Iterations........ 0 
C) Structure Loading: 
BX see bcedete hs edad ie aera rayne eres 2000.0000 
BY 5 ves diigo Si blew, Gb le ee eae wo STR ace 0.0000 
BM sete e ons tee woe Weare ae a leceie ose atcereae% 0.0000 
BA isso o eid eiadece sete aie to See eCS 0.0000 
D) Elemental Dimensions and Stress Distribution: 
Element Height Base Length 
1 2.00000 1.84550 8.00000 
2 2.00000 1.38407 8.00000 
3 2.00000 0.92261 8.00000 
4 2.00000 0.46154 8.00000 
E) Objective Function: 
Total structure Volume: 73.819420 
Node Stress 
1 §2018.20 
2 52020.44 
3 §2026.59 
4 52000.27 
5 9.4704286E-05 
F) Boundary Conditions: 
Node X-Displ Displ Slope 
1 1 1 1 
5 0 0 0 
G) Solution Vector: 
Node X-Displ Y-Displ Slope 
1 0.000000E+00 0.000000E+00 0.000000E+00 
2 0.508622E-01 0.221694E-08 -0.121376E-01 
: 3 0.197286E+00 0.860890E-08 -0.236977E-01 
4 0.433113E+00 0.189046E-07 -0.341030E-01 
5 0.742914E+00 0.324212E-07 -0.410363E-01 
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30000000.0 
52000.0 
4 


Volume 
29.52806 
22.14510 
14.76169 

7.38456 


OPTIMIZATION #1 


OPTIMIZATION SOLUTION 


ce ee ee a a a a a a a a a a ae a a a a a a a a en a A a a ww a aw a ww ww nr wr re ee ewe ere 


A) Problem Parameters: 


Arch Angle : 0.002 You 
Arch Radius: 1000000.000 Yie 
Arch Height: 2.000 No 


B) Derived Constants: 


No of System Nodal Points... 
No of Degrees of Freedom.... 
Length per Element.......... 


Number of Iterations 


C) Structure Loading: 


D) Elemental Dimensions and Stress Distribution: 


Element Height Base 
1 2.00000 1.84548 
2 2.00000 1.61572 
3 2.00000 1.38417 
4 2.00000 1.15351 
5 2.00000 0.92224 
6 2.00000 0.69289 
7 2.00000 0.46127 
8 2.00000 0.23665 
E) Objective Function: 
Total structure Volume: 66.495 
Node Stress 
1 §2015.55 
2 $1986.30 
3 52013.96 
4 52013.29 
5 52045.58 
6 $1955.51 
7 §2029.16 
8 §0706.71 
9 0.2237021 
F) Boundary Conditions: 
Node X-Displ Y-Displ 
1 1 1 
9 0 0 
G) Solution Vector: 
Node X-Displ Y-Displ 
1 0.000000E+00 0.000000E+00 
2 0.132929E-01 0.577893E-09 
3. 0.525036E-01 0.228824E-08 
4 0.117357E+00 0.511886E-08 
§ 0.207485E+00 0.905344E-08 
6 0.322357E+00 0.140683E-07 
7 0.461109E+00 0.201250E-07 
8 0Q.622199E+00 0.271538E-07 
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ngs Modulus: 
1d Strength: 
of Elements: 
9 
27 
4.0000 
1 
00.0000 
0.0000 
0.0000 
0.0000 
Length 
4.00000 
4.00000 
4.00000 
4.00000 
4.00000 
4.00000 
4.00000 
4.00000 
445 
Slope 
1 
0 
Slope 


0.000000E+00 
-0.650196E-02 
-0.129384E-01 
-0.192957E-01 
-0.255373E-01 
-0.316093E-01 
~0.373821E-01 
~0.425851E-01 


30000000.0 
52000.0 


Pwo) wo 


8 


Volume 
14. 
12. 
11. 
-22805 
-37795 
- 54310 
-69015 
. 89324 


76385 
92572 
07339 








9 0.801554E+00 0.349690E-07 -0.459655E-01 
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OPTIMIZATION #1 


OPTIMIZATION SOLUTION 


A) Problem Parameters: 


Arch Angle : 0.002 Youngs Modulus: 30000000.0 
Arch Radius: 1000000.000 Yield Strength: 52000.0 
Arch Height: 2.000 No of Elements: 12 
B) Derived Constants: 
No of System Nodal Points... 13 
No of Degrees of Freedom.... 39 
Length per Element..........- 2.6667 
Number of Iterations........ 4 
C) Structure Loading: 
EX cob hee bee oe We eo ON Sa es 2000.0000 
EY bos ceekiw acd ab Saas SES hee ece 0.0000 
BM pce into Ris ieee See we eS oS ore 0.0000 
BA 2 o6oSaince ecbee Be oo 8 suse bb weiss 0.0000 
D) Elemental Dimensions and Stress Distribution: 
Element Height Base Length Volume 
1 2.00000 1.84829 2.66667 9.85756 
2 2.00000 1.69537 2.66667 9.04196 
3 2.00000 1.54275 2.66667 8.22798 
4 2.00000 1.38887 2.66667 7.40732 
5 2.00000 1.23566 2.66667 6.59016 
6 2.00000 1.08215 2.66667 5.77148 
7 2.00000 0.93133 2.66667 4.96712 
8 2.00000 0.77510 2.66667 4.13389 
9 2.00000 0.63162 2.66667 3.36863 
10 2.00000 0.46759 2.66667 2.49383 
11 2.00000 0.31955 2.66667 1.70426 
12 2.00000 0.20000 2.66667 1.06667 
E) Objective Function: 
Total structure Volume: 64.630867 
Node Stress 
1 §1925.73 
2 51892.89 
3 51843 .36 
4 51329. 66 
5 51784.52 
6 $1739.70 
7 $1531.16 
8 51598.32 
9 50656.20 
10 51320.27 
11 50069.02 
12 39998.66 
13. 0.6707708 
F) Boundary Conditions: 
Node X-Displ Y-Displ Slope 
1 1 1 1 
13 0 0 0 


G) Solution Vector: 
Node X-Displ Y-Displ Slope 
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- 000000E+00 
-598324E-02 
-237427E-01 
-532194E-01 
-943462E-01 
-147043E+00 
-211205E+00 
-286684E+00 
.373299E+00 
-470718E+00 
-578546E+00 
-696051E+00 
.820672E+00 


gooooo0oo0oo0oo0o0°0o°c”c”e 


- O00000E+00 
-259434E-09 
-103343E-08 
- 231938E-08 
-411430E-08 
-641459E-08 
-921563E-08 
-125107E-07 
-162918E-07 
-205440E-07 
-252490E-07 
-303731E-07 
- 358010E-07 


=$1= 





0.000000E+00 
-0.442333E-02 
-0.882640E-02 
-0.132043E-01 
-0.175555E-01 
-0.218709E-01 
-0.261415E-01 
-0.303404E-01 
-0.344682E-01 
~0.384082E-01 
-0.422098E-01 
-0.455477E-01 
-0.473255E-01 





OPTIMIZATION #2 





OPTIMIZATION SOLUTION 


Arch Angle : 90.000 Youngs Modulus: 30000000.0 

Arch Radius: 32.000 Yield Strength: 52000.0 

Arch Height: 2.000 No of Elements: 12 
B) Derived Constants: 

No of System Nodal Points... 13 

No of Degrees of Freedom.... 39 

Length per Element.......... 4.1858 

Number of Iterations........ 1 


C) Structure Loadihg: 


BK 5 e oases Ss on Tories ne wife sevayal es ha save 0.0000 
PY seats Spend de 6 eins Rubin S ocw Wes SHES -2000.0000 * 
BM escdis sae esas ve. ae ei diw eit ayaa ose tehs 0.0000 
PA Gongs ne: oy 'b Fateat gfe Naps hS Over alison ey Doers 0.0000 
D) Elemental Dimensions and Stress Distribution: 
Element Height Base Length Volume 
1 2.00000 1.86405 4.18580 15.60506 
2 2.00000 1.84215 4.18580 15.42178 
3 2.00000 1.79464 4.18580 15.02403 
4 2.00000 1.71836 4.18580 14.38546 
5 2.00000 1.61213 4.18580 13.49611 
6 2.00000 1.47665 4.18580 12.36194 
7 2.00000 1.32535 4.18580 11.09533 
8 2.00000 1.16593 4.18580 9.76072 
9 2.00000 0.93161 4.18580 7.79905 
10 2.00000 0.71259 4.18580 5.96554 
11 2.00000 0.51683 4.18580 4.32666 
12 2.00000 0.36557 4.18580 3.06044 
E) Objective Function: 
Total structure Volume: 128.302124 
Node Stress 
1 519758.23 
2 51991. 84 
3 52001.67 
4 51952.61 
5 $1914.72 
6 51926.68 
7 §1573.75 
8 50482.76 : 
9 51888.91 
10 51924.89 
11 48451.66 « 
12 34653.70 
13 178.9231 
F) Boundary Conditions: 
Node X-Displ Y-Displ Slope 
1 1 1 1 
13 0 0 0 
G) Solution Vector: 
Node X-Displ Y-Displ Slope 
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qoooo0oo0ooo0o°co0o0o 


. 000000E+00 
-149422E-01 
-589253E-01 
-129477E+00 
-222644E+00 
.333172E+60 
.454693E+00 
.579860E+00 
. 700496E+00 
- 808353E+00 
.895118E+00 
-952352E+00 
.972583E+00 


- 000000E+00 
-105421E-02 
-987875E-02 
.339057E-01 
.799319E-01 
.153870E+00 
- 260536E+00 
-403368E+00 
-584032E+00 
-802894E+00 
-105869E+01 
.134670E+01 
-165575E+01 
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0. 
-0. 
-0. 
-0. 
-0. 
-0. 
-0. 
-0. 
-0. 
~0. 
~0. 
~0. 
-0. 


000000E+00 
714707E-02 
142563E-01 
213028E-01 
282739E-01 
351633E-01 
419651E-01 
486117E-01 
549785E-01 
613227E-01 
673512E-01 
723963E-01 
747875E-01 














OPTIMIZATION #3 


OPTIMIZATION SOLUTION 


A) Problem Parameters: 


Arch Angle : 90.000 Youngs Modulus: 30000000.0 

Arch Radius: 32.000 Yield Strength: 52000.0 

Arch Height: 2.000 No of Elements: 12 
B) Derived Constants: 

No of System Nodal Points... 13 

No of Degrees of Freedom.... 39 

Length per Element.......... 4.1858 

Number of Iterations........ 1 


BX ives es Sais eee Sieiee aw 8 oe So Nie 2000.0000 
BY. een Sk Se Baw bac ated eee oh 0.0000 
BM yo ace tecave diei toc wieisibie biel Gum Sia-eceee 0.0000 
PA iis eases Seb. 8-8 beers, os EMSs ae 0.0000 
D) Elemental Dimensions and Stress Distribution: 
Element Height Base Length Volume 
1 2.00000 1.84614 4.18580 15.45518 
2 2.00000 1.68690 4.18580 14.12203 
3 2.00000 1.43532 4.18580 12.01594 
4 2.00000 1.19577 4.18580 10.01052 
5 2.00000 0.97029 4.18580 8.12291 
6 2.00000 0.75740 4.18580 6.34069 
7 2.00000 0.56761 4.18580 4.75181 
8 2.00000 0.40093 4.18580 3.35641 
9 2.00000 0.26110 4.18580 2.18585 
10 2.00000 0.20000 4.18580 1.67432 
11 2.00000 0.20000 4.18580 1.67432 
12 2.00000 0.20000 4.18580 1.67432 
E) Objective Function: 
Total structure Volume: 81.384300 
Node Stress 
1 52001.90 
2 51862.35 
3 §1955.02 
4 51942.95 
5 $1853.82 
6 51974.19 
7 §1925.02 
8 51871.48 
9 51649.30 
10 41477.37 
11 21322.36 
12 9098.637 
13 4994.856 
F) Boundary Conditions: 
Node X-Displ Y-Displ Slope 
1 1 1 1 
13 0 0 0 


G) Solution Vector: 
Node X-Displ Y-Displ Slope 
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0.000000E+00 
0.144840E-01 
0.557774E-01 
0.120917E+00 
0.206307E+00 
0.307129E+00 
0.417585E+00 
0.531063E+00 
0.640343E+00 
0.737912E+00 
0.813769E+00 
0.860985E+00 
0.877143E+00 


0.000000E+00 
-0.944371E-03 
-0.914169E-02 
-0.312207E-01 
-0.732726E-01 
-0.140544E+00 
-0.237250E+00 
-0.366366E+00 
-0.529394E+00 
-0.726163E+00 
-0.948447E+00 
-0.118418E+01 
-0.142556E+01 


ee 


0.000000E+00 
-0.677752E-02 
-0.131680E-01 
-0.195029E-01 
-0.257570E-01 
-0.319051E-01 
-0.379512E-01 
-0.438425E-01 
-0.495298E-01 
-0.549162E-01 
-0.570368E-01 
-0.577029E-01 
-0.578414E-01 





OPTIMIZATION #4 


OPTIMIZATION SOLUTION 


ie we ne we we we ee ew wm — a a a oo ww a ww ww ww ww we wr we re re ower 


A) Problem Parameters: 





Arch Angle : 90.000 Youngs Modulus: 30000000.0 
Arch Radius: 32.000 Yield Strength: §2000.0 
Arch Height: 2.000 No of Elements: 12 
B) Derived Constants: 
No of System Nodal Points... 13 
No of Degrees of Freedom.... 39 
Length per Element.......... 4.1858 
Number of Iterations........ 1 : 
C) Structure Loading: 
MK Seis: ha Wiener erbuel er ei eeerd do eae BRS 0.0000 
BY oie ca cen eae tees 0.0000 : 
BM iis 655 Save eietels tree Bis eels = a0 0.0000 
BA Sire ecb Me ea heise eee Le hiere -100.0000 
D) Elemental Dimensions and Stress Distribution: 
Element Height Base Length Volume 
1 2.00000 2.08018 4.18580 17.41440 
2 2.00000 1.93027 4.18580 16.15947 
3 2.00000 1.76172 4.18580 14.74842 
4 2.00000 1.55158 4.18580 12.98917 
5 2.00000 1.31450 4.18580 11.00444 
6 2.00000 1.06716 4.18580 8.93384 
7 2.00000 0.82127 4.18580 6.87536 : 
8 2.00000 0.59123 4.18580 4.94951 
9 2.00000 0.39007 4.18580 3.26550 
10 2.00000 0.22180 4.18580 1.85686 
11 2.00000 0.20090 4.18580 1.67432 
12 2.00000 0.20000 4.18580 1.67432 
E) Objective Function: 
Total structure Volume: 101.545610 
Node Stress 
1 $2103.21 
2 52039.99 
3 §2033.59 
4 §2017.27 
5 §2015.93 
6 51958.49 
7 §1950.00 
8 51920.60 4 
9 51610.81 
10 §1626.45 
11 25188. 80 . 
12 §614.041 
13 1096.744 
F) Boundary Conditions: 
Node X-Displ Y-Displ Slope 
1 1 1 1 
13 0 0 0 
G) Solution Vector: 
Node X-Displ Y-Displ Slope 








OBDIARMN &®WNHH 


0.000000E+00 
0.145714E-01 
0.573794E-01 
-125999E+00 
-216392E+00 
- 323304E+00 
- 440430E+00 
- 560627E+00 
-676171E+00 
-778995E+00 
-861090E+00 
- 913642E+00 
- 931478E+00 


oooooocncoooo 


- 000000E+00 
.111348E-02 
-978175E-02 
.332242E-01 
- 779492E-01 
-149537E+00 
- 25241 0E+00 
- 389635E+00 
-562738E+00 
-771440E+00 
-101350E+01 
.127777E+01 
-154990E+01 
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0.000000E+00 
-0.69174.E-02 
-0.138459E-01 
-0.206572E-01 
-0.273459E-01 
-0.339032E-01 
-0.403032E-01 
-0.465215E-01 
-0.525061E-01 
-0.581423E-01 
-0.633131E-01 
-0.650579E-01 
-0.651839E-01 














OPTIMIZATION #5 


OPTIMIZATION SOLUTION 


cca ewan cee me eee ee awe eee ee a Se ee 6 eS 8 8S 8 0 


A) Problem Parameters: 


Arch Angle : 90.000 Youngs Modulus: 
Arch Radius: 32.000 Yield Strength: 
Arch Height: 2.000 No of Elements: 
B) Derived’ Constants: 
No of System Nodal Points... 13 
No of Degrees of Freedom.... 39 
Length per Element..........- 4.1858 
Number of Iterations........ 
C) Structure Loading: 
EX cos SAS is ow ese ole aioe 0.0000 
BY. osc bebe erbiie ws eia © bi wiseye sieeve -8000.0000 
BM oo acece oats 2 Sse ero oreue ovate ace. este 0.0000 
BA bciice's oo.ccdle ele areal erede 0 eae 8 0.0000 
D) Elemental Dimensions and Stress Distribution: 
Element Height Base Length 
1 2.00000 0.51296 4.18580 
2 2.00000 0.90694 4.18580 
3 2.00000 1.11658 4.18580 
4 2.00000 1.47904 4.18580 
5 2.00000 1.19599 4.18580 
6 2.00000 1.12703 4.18580 
7 2.00000 1.53118 4.18580 
8 2.00000 0.57650 4.18580 
9 2.00000 0.59349 4.18580 
10 2.00000 1.34594 4.18580 
11 2.00000 2.20595 4.18580 
12 2.00000 3.12364 4.18580 
E) Objective Function: 
Total structure Volume: 131.561722 
Node Stress 
1 8078.639 
2 52000.93 
3 50895.98 
4 52005. 46 
5 52048.33 
6 §2010.10 
7 42175.72 
8 49131.36 
9 4933.140 
10 52013.68 
11 §2043.57 
12 §2041.75 
13 52027.50 
F) Boundary Conditions: 
Node X-Displ Y-Displ Slope 
1 1 1 0 
13 1 0 1 
G) Solution Vector: 
Node X-Displ Y-Displ Slope 
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30000000.0 
52000.0 


12 


Volume 
- 29430 
-59253 
.34755 
- 38193 
.01239 
- 43503 
. 81847 
- 82624 
- 96845 
11.26764 
18.46734 
26.14985 


we pa} 
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0.000000E+00 
-0.823069E-01 
-0.143411E+00 
-0.178618E+00 
-0.189389E+00 
-0.177714E+00 
~0.146174E+00 
-0.105000E+00 
-0.619691E-01 
-0.273824E-01 
-0.803875E-02 
-0.764124E-03 

0.000000E+00 


'#petetepetud 
oooooocoocooo0°0o 


. 000000E+00 
-426506E-02 
-157328E-01 
-270851E-01 
.319111E-01 
-234612E-01 
-494672E-02 
-525035E-01 
.118718E+00 
-190906E+00 
-249018E+00 
-286585E+00 
- 300012E+00 
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0.208648E-01 
0.173341E-01 
0.118840E-01 
0.554722E-02 
0.244944E-04 
-0.683719E-02 
-0.132138E-01 
-0.165616E-01 
-0.201392E-01 
~0.168476E-01 
-0.117570E-01 
~0.606703E-02 
0.000000E+00 








OPTIMIZATION #6 


OPTIMIZATION SOLUTION 


we ee ee ea ow a ee a a a a a a ow we a a a ww ww ww we wr wn ww wr ww www wr wre wre 


Arch Angle : 90.000 Youngs Modulus: 30000000.0 

Arch Radius: 32.000 Yield Strength: 52000.0 

Arch Height: 2.000 No of Elements: 12 
B) Derived Constants: 

No of System Nodal Points... 13 

No of Degrees of Freedom.... 39 

Length per Element.......... 4.1858 

Number of Iterations........ 1 


C) Structure Loading: 


1D Cee i eT ee ee eee ae ee 0.0000 
BY. ie Fie eo Riw aire aisha eee as Wie whee ues -8000.0000 
BM ooo eeterghate side leliel eels w: atbceeecdvone 0.0000 
FA o:dioo S056 Ss oi0 Sie 0: ose ha ete ce-e evans 0.0000 
D) Elemental Dimensions and Stress Distribution: 

Element Height Base Length Volume 

1 2.00000 1.72707 4.18580 14.45838 

2 2.00000 0.94577 4.18580 7.91762 

3 2.00000 0.27194 4.18580 2.27661 

4 2.00000 0.65040 4.18580 5.44488 

5 2.00000 1.09075 4.18580 9.13130 

6 2.00000 0.78986 4.18580 6.61236 

7 2.00000 1.06223 4.18580 8.89255 

8 2.00000 0.75388 4.18580 6.31116 

9 2.00000 0.25849 4.18580 2.16400 

10 2.00000 0.93275 4.18580 7.80861 

11 2.00900 1.75263 4.18580 14.67236 

12 2.00000 2.65143 4.18580 22.19675 


E) Objective Function: 
Total structure Volume: 107.886574 
Node Stress 


51997.34 
51279.07 
11 $2434.15 
12 51908.94 
13 51935.26 


F) Boundary Conditions: 


Node X-Displ Y-Displ Slope 
1 1 1 1 
13 1 0 1 


G) Solution Vector: 
Node X-Displ Y-Displ Slope 
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0.000000E+00 
-0.121854E-01 
~0.442765E-01 
-0.861533E-01 
-0.112986E+00 
-0.121058E+00 
-0.110454E+00 
-0.850803E-01 
-0.537900E-01 
-0.252856E-01 
~0.744947E-02 
-0.654979E-03 

0.000000E+00 


0.000000E+00 
0.457452E-03 
0.615066E-02 
0.177060E-01 
0.297132E-01 
0.342979E-01 
0.237526E-01 
-0.622361E-02 
-0.547404E-01 
-~0.118485E+00 
-0.173190E+00 
-0.209063E+00 
~0.222069E+00 
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0.000000E+00 
0.527930E-02 
0.9475S56E-02 
0.907372E-02 
0.435047E-02 
~0.109393E-03 
-0.701848E-02 
-0.115516E-01 
-0.155000E-01 
-0.155501E-01 
~-0.110976E-01 
-0.580765E-02 
0.000000E+00 
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APPENDIX E 


ARCH_OPT.FOR SOURCE CODE 


PROGRAM ARCH OPTIMIZATION 
* 


ARCH OPTIMIZATION ANALYSIS CODE 7 
* 





ALPHA....TRANSFORMATION ANGLE OF ELEMENT (ANGLE TO X~-AXIS) 

ANGLE....TOTAL ANGLE OF ARCH (IN DEGREES) 

BASE..... DOT ARRAY CONTAINING THE ELEMENTAL BASE DIMENSIONS 

BASEL....DOT ARRAY CONTAINING THE ELEMENTAL BASE DIMENSIONS LOWER 
SIDE CONSTRAINT 

BASEU....DOT ARRAY CONTAINING THE ELEMENTAL BASE DIMENSIONS UPPER 
SIDE CONSTRAINT 


BETA ....TRANSFORMATION ANGLE OF ELEMENT (ANGLE TO Y-AXIS) 
B_l...... BOUNDARY TERMS APPLIED AT END "1" 
BQ. sees BOUNDARY TERMS APPLIED AT END "2" 
Cl,..,C5...CONSTANTS RELATED TO ELEMENT STIFFNESS COEFFICIENTS 
CLAN..... CONCENTRATED LOAD APPLICATION NODE (THE NODE FX,FY,FM ARE 
APPLIED 
COUNT....COUNTS THE NUMBER OF ITERATIONS COMPLETED 
DOF......DEGREE OF FREEDOMS (UNKNOWN DISPLACEMENTS & SLOPES) ‘ 
DSN...... DESIGN VARIABLE FOR EACH ELEMENT , 


DV1BG....DESIGN VARIABLE #1 (BASE DIMENSION) INITIAL ESTIMATE 
DV1LO....DESIGN VARIABLE #1 (BASE DIMENSION) LOWER SIDE CONSTRAINT 
DV1UP....DESIGN VARIABLE #1 (BASE DIMENSION) UPPER SIDE CONSTRAINT 


EK......-. 6X6 ELEMENT STIFFNESS MATRIX IN LOCAL X,Y COORDINATES 

EKPR..... 6X6 ELEMENT STIFFNESS MATRIX IN ELEMENT LOCAL COORDINATES 

ELEN..... LENGTH OF ELEMENT 

| as FORCE VECTOR OF SYSTEM 

FAL. wees CONSTANT DISTRIBUTED LOAD IN X DIRECTION FROM END TO END 

FM. ...64- CONCENTRATED MOMENT AT FREE END 

FRX. wee eee CONCENTRATED LOAD IN X DIRECTION AT FREE END 

FY. ....6. CONCENTRATED LOAD IN Y DIRECTION AT FREE END 

Giwsecaes THE ARRAY OF CONSTRAINT FUNCTIONS 

GAMMA....6X6 ELEMENT TRANSFORMATION MATRIX 

GK....... (NDOF)X(NDOF) GLOBAL STIFFNESS MATRIX 

Heo wseveces DEPTH OF ARCH SECTION 

INDSN....INITIAL (UNIFORM) DESIGN DIMENSION se 
INFO..... DOT PARAMETER USED TO SIGNAL THAT THE OPT IS COMPLETE 

IPRINT...DOT PARAMETER USED SELECT THE DATA OUTPUT FORMAT 

IPRM..... DOT SELECTABLE INTEGER PARAMETERS : ° 


ITERATE..THE NUMBER OF TIMES DOT IS TO BE RELOADED WITH THE 
PRECEEDING DATA 

IWK...... DOT INTERNAL WORK SPACE ARRAY 

METHOD...DOT PARAMETER USED TO DEFINE THE OPTIMIZATION METHOD 

MINMAX...DOT PARAMETER USED TO MINIMIZE/MAXIMIZE THE PROBLEM 


NCON..... NUMBER OF DESIGN CONSTRAINTS 

NDOF..... NUMBER OF DEGREES OF FREEDOM 

NEL... 0.5.3 NUMBER OF ELEMENTS 

NRIWK....DOT INTERNAL WORK SPACE ARRAY DIMENSION 

NRWK..... DOT INTERNAL WORK SPACE ARRAY DIMENSION 
-~92- 











NSNP.....NUMBER OF SYSTEM NODAL POINTS 

OBJ......THE OBJECTIVE FUNCTION OF THE OPTIMIZATION 
OPTDCS...OPTIMIZATION DECISION TO OPTIMIZE THE PROBLEM OR NOT 
P1...P4..PARAMETER DIMENSION CORRESPONDING TO THE NEL, NSNP, NCON, 
AND NDOF, RESPECTIVELY 

PHT ices SUBTENDED ELELENT ANGLE (ALSO, PHIANG IN DEGREES) 
PRCSN....THE PRECISION DESIRED TO SOLVE THE FEM SYSTEM OF EQUATIONS 
RADIUS...ARCH RADIUS 

RPRM..... DOT SELECTABLE REAL PARAMETERS 

SIGMA_B..THE ELEMENTAL NORMAL STRESS DUE TO BENDING 

SIGMA_N..THE ELEMENTAL NORMAL STRESS DUE TO AXIAL FORCES 
SIGMA_T..THE MAXIMUM TOTAL STRESS IN EACH ELEMENT 

Sie eueres ..THE "DISPLACEMENT" VECTOR OF THE SYSTEM OF LINEAR EQUATIONS 
WK Secere DOT INTERNAL WORK AREA 

X.......-GLOBAL HORIZONTAL COORDINATE 

Me 2d.d ie ayeece GLOBAL VERTICAL COORDINATE 

YIELD....YIELD STRENGTH OF THE ARCH MATERIAL 

YOUNG. ...YOUNG’S MODULUS OF THE ARCH MATERIAL 

eK KKK KKK RE KKKAKKEKKEKKKEEKE EERE KEKKKKKKAKKKKKEKKEKKKKKKKEKEK 
....@eclare the variables............. Base) sesso we heete ei te daves'e.© eee 
INCLUDE ‘ ARCH_COM. FOR’ 


AA Aweeene een ee_ee ee eee ERE 
GQ 


~...read the input parameters... .... cee cee eee eee eee e nes 
OPEN (8, FILE=’ARCH_IN.DAT’, STATUS=’OLD’ ) 
READ (8,*) ANGLE, RADIUS, YOUNG, YIELD, NEL, METHOD, IPRINT, DV1BG, 


& DV1LO, DV1UP,H, CLAN, FX, FY, FM, FA, OPTDCS, ITERATE, PRCSN, 
& BX1, BY1, BM1, BX2, BY2, BM2, LABEL 
c 
Cc ....define constants.......... sna duo dua eicaraterepa ee 5 eiiag: 8 fevsc evga tere vars. 2 
NCON = 3*NEL 
NSNP = NEL + 1 
NDOF = 3*NSNP 
Cc 
Cc ....determine the system nodal coord and element orientation.. 
CALL GEOMETRY (NEL, NSNP, ANGLE, RADIUS, X,Y, ALPHA, BETA, ELEN) 
c 
Cc ....define the size of the work arrays for DOT............... F 
NRWK = 38800 
NRIWK = 400 
Cc 
c ....optimize the problem..... steer arels cvleceershe de tos eer tarpatsarimaet oles Nets fe\'ey'sh 6 
CALL OPTIMIZATION_TOOL 
c 
Cc -...compile and format the output.......... 0. ce eee wenes sia be cata 
CALL ARCH_OUTPUT 
Cc 
END 
PERE RPE SERSRESSRESZASRSSSSSESALESELSESSSESESASESSESE LESSEE SES SESE SEER SSS ESS SSS 
* * 
SUBROUTINE GEOMETRY (NEL, NSNP,ANGLE, RADIUS, X, ¥, ALPHA, BETA, ELEN) 
Cc ES ee EE ETT SRST ETC SEES SLES AAATTETSETTCETCTCTETETC ESAS KBESEC ESS SS 
c | This routine is used by main ARCH_OPTIMIZATION to generate | 
c | the x-, y-coordinates of each system node, to determine | 
Cc | the orientation of each element, andto calculate the | 
c | length of each element. | 
Cc oe Se Sr er St ee ee ee 2 ee Sea ee 
Cc ....declare the variables. ....... ccc ccc cc cee eee ee ce teens 


INTEGER NEL, NSNP,P1,P2 

PARAMETER (P1#=#32,P2=33) 

REAL ANGLE, RADIUS, ELEN, X(P2),Y¥(P2),ALPHA(P1) BETA(P1), 
é PI, PHI,ANG, YNUM, XDEN 
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aa 


100 





PARAMETER(PI = 3.141593) 


....determine the geometric constants... .... ccc e cence ene rnes 
PHI = (ANGLE/NEL) * (PI/189.0) 


X(1) = 0.0 
¥(1) = 0.0 


DO 100 i=2, NSNP 

ANG = (i-1.0)*PHI 

X(i) = RADIUS * (1.0 - COS (ANG) ) 

Y(i) = RADIUS * SIN(ANG) 

YNUM = (Y(i) - Y¥(i-1)) 

XDEN = (X(i) - X(i-1)) 

ALPHA(i-1) = ATAN2 (YNUM, XDEN) 

BETA(i-1) = (PI/2.0) - ALPHA(i-1) 

CONTINUE 


...-determine the length of each element...................--- 
ELEN = SORT (X(2)**2.0 + Y(2)**2.0) 


RETURN 
END 


RH KR RE KHMER KKK KKK KKK REE KR EKEKKEKKKEKKKKKKEKKEKKKKKKEK 
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aqaaQq aan aa aaAgNgNAN 


aq an 


aa 


100 





* 
SUBROUTINE OPTIMIZATION_TOOL 
This subroutine directs the program flow optimization decision 
i.e., optimize the problem or not. It also serves to set up & 
execute the DOT optimization software. 
SETRCCCCCCCKCCSCKLELECKC CECT SESSA CSTE eT SSE SESSA SASS See 
..-..declare the variables.......... ia alaetacesaneeoee er ae a are ore rah 
INCLUDE ’ARCH_COM.FOR’ 
INTEGER i 


-+..zero out the RPRM and IPRM arrayS.......... cece cece ec enee 
DO 100 i1,20 

RPRM(i) = 0.0 

IPRM(i) = 0 
CONTINUE 


oes PNLELALS Ze" COUNT eee a ois bi gite valor al oS ore: Berea Se ew Seve nelee eh e:.6 Whe yee 15) eae 
COUNT = 0 
....refine the constraint tolerence......... ccc cece eee ween 


RPRM(2) = 0.0001 
RPRM(3) = 0.0001 


-...turn off DOT’s auto scaling... .... ccc cee ewe cece reece cen eee 
IPRM(2) = -1 


-...increase DOT’s default number of iterations............... 
IPRM(3) = 1000 
IPRM(8) = 1000 


....increase DOT’s number of consecutive convergence criteria. 
IPRM(4) = 3 
IPRM(9) = 3 


....define MINMAX=-1 to minimize the objective function....... 
MINMAX = -1 

















Cc 
Cc ....dimitialize the design variable limits and best guess...... 
DO 200 i=1,NEL 
BASE (i) = DV1BG 
BASEL(i) = DV1LO 
BASEU(i) = DV1UP 
200 CONTINUE 
Cc 
Cc -..-Make optimization decision... ..... cee cee eee eee eens 
IF (OPTDCS .NE. 1) THEN 
CALL EVAL 
RETURN 
ENDIF 
Cc 
P Cc es. keaGy tO Optimize... 2s. sce o ee we eee we ces Ba caitetwis Sieyloe- Sa cerec eee 
INFO = 0 
Cc 
° 300 CALL DOT (INFO, METHOD, IPRINT, NEL, NCON, BASE, BASEL, BASEU, OBJ, 
& MINMAX, G, RPRM, IPRM, WK, NRWK, IWK, NRIWK) 
Cc 
Cc -..-evaluate the objective function and constraints........... 
IF (INFO .GT. 0) THEN 
CALL EVAL 
GOTO 300 
ENDIF 
Cc 
Cc ....refine the solution vector by reoptimizing................ 
IF (COUNT .LT. ITERATE) THEN 
INFO = 0 
COUNT = COUNT+1 
GOTO 300 
ENDIF 
c 
RETURN 
END 


(SERS REREASELASASRERARLASESLELSESSAESESESLESSESESSESESELESE SEAS SESS SSAS SEES SES 
* x 


SUBROUTINE EVAL 


Cc (St et A Oe OG a a ee ES eS SS 

Cc This subroutine is used to evaluate the Objective function, 

Cc constraint functions, and side constraints of the optimization 

Cc problem. 

(es ae Se ee ee Se See ee SS Be See a 

Cc ..-.declare the variables... ..... ccc ccc ce eee ee eens 
INCLUDE ‘ARCH_COM.FOR’ 
INTEGER i,j 

Cc 

as Cc ....calculate the objective function.................26- pea ates 

OBJ = 0.0 

Cc 

" DO 100 i=1,NEL 
OBJ = OBJ + BASE(i)*H*ELEN 
100 CONTINUE 

Cc 

Cc ....initialize the design constraint vector............. sale, anes 
DO 200 i=1,NCON 

G(i) = 0.0 
200 CONTINUE 
Cc 
Cc ....determine the design constraints......... 2... cece ee ee ee eee 


CALL ARCH_STRESS 
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c 
DO 210 i=1,NEL 
IF (SIGMA T(i) .GE. SIGMA_T(i+1)) THEN 
SIGMA = SIGMA_T (i) 
ELSE 
SIGMA = SIGMA _T(it1) 
ENDIF 
c 


G(i) = (SIGMA/YIELD) - 1.0 
210 CONTINUE 


DO 220 i*1,:° EL 
}=i+NEL 
G(j) = (BASE(i)/(3.0*H)) - 1.0 
220 CONTINUE 


DO 230 i=1,NEL 
#i+(2*NEL) 
G(j) = H/(10.0*BASE(i)) - 1.0 ry 
230 CONTINUE 


RETURN 
END 


RK KKK KKK KE KKK KK KKK KKK KK KEK KEK RE KK KKK KK KE KEKE KKEEKEKEKEEKKKKEKK 
* * 
SUBROUTINE ARCH STRESS 

‘his subroutine is used to perform the Finite Element analysis 
of the stresses developed in an arch or beam for a given load- 
ing. 

BERR RSSSETAECKLESRSEALLoOS SERBS OCS SS SS SSS SSS SSS SSS SES SSS SE ETT 
--.-declare the variables.......... cc eee eee eee syela. sie neie se 
INCLUDE ’ARCH_COM.FOR’ 

INTEGER IPVT (99) 

REAL GK (P4,P4),F (P4) 

REAL*8 BK(P4,P4),BF (P4),BU(P4) ,FAC(9801) , WORK (99) 


aaqaqaa 


-...form the eiement and system matrices.......... cece cece eee 
CALL FORM (NEL, NDOF, ALPHA, BETA, H, ELEN, YOUNG, BASE, GK) 


-...form the Force vector, F.... ccc cee wc wee nce tee eee wees 
CALL FORCE_VECTOR (NEL, NDOF, ELEN, ALPHA, BETA, FA, F) 


-...set the boundary conditiona and loads................20020: 
CALL BNDARY (NDOF, GK, CLAN, FX, FY, FM, F,BX1,BY1,BM1, BX2, BY2, BM2) 


-...SOlve the system 01 equations........... ce ee ec eee eee . 
IF (PRCSN .EQ. 2) THEN 
...-Change GK and F arrays to double precision........ 
CALL UPSCALE (NDOF,GK,F, BK, BF) e 
...-SOlve the system of equationsS............ cece eeeee 
CALL DLZARG (NDOF, BK, P4, BF,1,BU, FAC, IPVT, WORK) 
....Cchange BU array to single presicion............... 
CALL DOWNSCALE (NDOF, BU, U) 
ELSE 
Cc ....SOlve the system of equations............. cece wees 
CALL L2ARG (NDOF,GK,P4,F,1,U, FAC, IPVT, WORK) 
ENDIF 


oO aq a ag AN AN AA 


Cc ....determine the stress distribution........... 0. ccc eee ee wees 
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CALL STRESS (X,Y, ALPHA, BETA, U, NEL, ELEN, YOUNG, H, SIGMA _ T) 
Cc 

RETURN 

END 


KK KKK KH KKK KKK KEK KKK KEE REREKKKKEREKRREKKEKEKEEKKEKKKKKKEKEKKKKKKKKK 

* «x 
SUBROUTINE FORM (NEL, NDOF, ALPHA, BETA, H, ELEN, YOUNG, BASE, GK) 

Cc DEE SE A Se Se Se AD 

Cc This subroutine is used to construct the global stiffness mat- 

Cc rix for the arch problem. 

Cc ot fe ee ee 2 de See 2 ee Ee 8 eee ee 8 ee ee eS ee ee ee 

Cc sees Geclare the variables soi oeies ieee sek ele os be cele eeee ele ee ae we’ 
INTEGER NEL, NDOF,NPOW, IEL,1I,J,K,II,JJ,KK, III, JJJ,P1,P2,P4 

Cc 

Cc 


PARAMETER (P1=32, P2=33, P4=99) 


REAL ELEN, H, BASE (P1) ,ALPHA(P1), BETA(P1) , YOUNG, 
C1,C2,C3,C4,C5, CA, CB, EKPR(6, 6) ,GAM(6,6),EK(P1,6,6), 
GAMMA (P1,6,6),GK(P4,P4) ,EKGA(6, 6) , GAEKGA (6, 6), 
ALPHAI, BETAI 


HM wm 


agQ 


....define the constantsS Cx... .. ccc ccc cee eee ee cere ee ee eeee 
NPOW = 1 

Cl = YOUNG*H/ELEN 

C2 = (H/ELEN) **2.0 

C3 = (H**2.0)/(2.0*ELEN) 

C4 = (H**2.0)/3.0 

cS = C4/2.0 


Cc ....initialize the EKPR and GAM arrayS..............eeeeeeeees 
DO 100 II = 1,6 
DO 90 JJ= 1,6 
EXPR(II,JJ) = 0.U 
GAM(II,JJ) = 0.0 
90 CONTINUE 
100 CONTINUE 


aa 


....initialize the EK and GAMMA arrayS..............cecceeeees 
DO 130 IEL = 1,NEL 
DO 120 I = 1,6 
po 110 J = 1,6 
EK(IEL,i,J) = 0.0 
GAMMA(IEL,I,J) = 0.0 
110 CONTINUE 
120 CONTINUE 
130 CONTINUE 


ana 


....determine the EKPR matrix....... cece ec ecw ee eee ee eee eeeee 
EKPR(1,1) c1 

EKPR (1, 4) 
EKPR (2,2) 
EKPR (2, 3) 
EKPR (2,5) 
EKPR (2, 6) 
EKPR (3, 2) 
EKPR (3, 3) 
EKPR (3,5) 
EKPR (3, 6) 
EKPR (4,1) 
EKPR (4, 4) 
EKPR (5, 2) 


Oeeeeannee @ ban 
Q 
~ 
» 
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EKPR(5,3) = -Ci*C3 
EKPR(5,5) = C1*C2 
EKPR(5,6) = -C1*C3 
EKPR(6,2) = C1*C3 


EKPR(6,3) = C1*C5 
EKPR(6,5) = -C1*C3 
EKPR(6,6) = C1*C4 


....initialize the GK array... .... cc ccc cee ccc eee tert ewe eres 
DO 150 I = 1, NDOF 
DO 140 J = 1, NDOF 
GK(I,J) = 0.0 
140 CONTINUE 
150 CONTINUE 


....Getermine the GAMMA matrix... .. ccc ccc eee cece tcc eens 
DO 170 IEL = 1,NEL 

ALPHAI = ALPHA(IEL) 

BETAI = BETA(IEL) 

CA = COS (ALPHAT) 

CB = COS (BETAI) 


GAMMA(IEL,1,1) = CA 
GAMMA(IEL,1,2) = CB 
GAMMA (IEL,2,1) = -CB 
GAMMA (IEL,2,2) = CA 
GAMMA (IEL,3,3) = 1.0 
GAMMA (IEL,4,4) = CA 
GAMMA(IEL,4,5) = CB 
GAMMA (IEL,5,4) = -CB 
GAMMA(IEL,5,5) = 


¢ 


GAMMA (IEL, 6, 6) 
170 CONTINUE 


....initialize the EKGA and GAEKGA arrays....... epee evecare hie 3% 
DO 270 IEL = 1, NEL 
DO 190 III = 1,6 
DO 180 JJJ = 1,6 
EKGA(III,JJJ) = 0.0 
GAEKGA(III,JJJ) = 0.0 
180 CONTINUE 
190 CONTINUE 


...determine the EKGA array.......... ccc eee cece cece ec erences 
DO 220 I = 1,6 
DO 215 J = 1,6 
DO 210 K = 1,6 
EKGA(I,J) = EKGA(I,J) + EKPR(I,K) *GAMMA(IEL,K, J) 

210 CONTINUE 
215 CONTINUE 
220 CONTINUE 


....determine the GAEKGA array....... 2 cece ccc n cece ence cr eceee 
DO 235 J = 1,6 
DO 230 K = 1,6 
GAEKGA (I,J) = GAEKGA(I,J) 


& +GAMMA (IEL, K, I) *EKGA (K, J) 
230 CONTINUE 
235 CONTINUE 
240 CONTINUE 
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Cc ....copy the GAEKGA array into the EK array.......-----seseee- 
DO 260 I = 1.6 
DO 250 J = 1,6 
EK(IEL,I,J) = GAEKGA(I,J) 
250 CONTINUE 
260 CONTINUE 
270 CONTINUE 


Cc ~.--construct the GK matrix... ccc cc ee ew ee eee eee eee eee 
DO 300 IEL = 1, NEL 
II = 3* (IEL-1) 
Do 290 J = 1, 6 
JJ=II+ J 
DO 280 K = 1, 6 
KK = II + K 
GK(JJ,KK) = GK(JJ,KK) 
& +EK (IEL, J, K) * (BASE (IEL) * *NPOW) 
280 CONTINUE 
290 CONTINUE 
300 CONTINUE 


RETURN 
END 


PK MIM HHH HK KK HII KKK KKK KKK KKK KKK KKK KEKE KKK KEKKKEKKEKKKKKKK 
* « 
SUBROUTINE FORCE_VECTOR (NEL, NDOF, ELEN, ALPHA, BETA, FA, F) 

REECE SKSRESE SACs TESESE ALBEE SSS TeEE CTI sEeEseTT 
This subroutine is used to construct the force vector for the 
FEM problem specified. 


Rr ate ee A eae A Sate a Te eS eS aes 


.... declare the variables.......... eetaseres jibe Sueseietecere Wise a.sidwlere 
INTEGER NEL, NDOF,i,1I1,12,13,P1,P4 


PARAMETER (P1=32, P4=99) 
REAL ELEN, ALPHA (P1),BETA(P1),FA,F(P4) 


aa a a anaana 


so ee LOEM: CHE: F—VeECC OR 6c ces. 8s ea wkd wea 6 Saree elie! Wis. 6 088 we SES Be wie oe ae 
F(1) = (ELEN/2.0) * (-COS(BETA(1))) 

F(2) = (ELEN/2.0) * (COS (ALPHA(1))) 

F(3) = (ELEN**2.0)/12.0 


DO 100 i=2,NEL 
Il = (i-1)*3 +1 
I2 = (i-1)*3 + 2 
I3 = (i-1)*3 + 3 


F(I1) = (ELEN/2.0) * (-COS (BETA (NEL) ) ) 
& + (ELEN/2.0) * (-COS (BETA (NEL-1) ) ) 
F(I2) = (ELEN/2.0) * (COS (ALPHA (NEL) ) ) 
& + (ELEN/2.0) * (COS (ALPHA (NEL-1) ) ) 
F(I3) = 0.0 
100 CONTINUE 


F(NDOF-2) = (ELEN/2.0) * (-COS (BETA (NEL) ) ) 
F(NDOF-1) = (ELEN/2.0) * (COS (ALPHA (NEL) ) ) 
F (NDOF) = -(ELEN**2.0)/12.0 


Cc ~...scale the F-vector by FA... .. 2. ccc ccc cee cece ener rene ee reee 
DO 200 i=1,NDOF 
F(i) = FA*F (i) 
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200 CONTINUE 
Cc 
RETURN 
END 


WHEE KKK KKH KEKREKEKEKKEKEKEKKKK KKK KEKEKKEKEKKERHEEEEKKKKKKK 
* * 
SUBROUTINE BNDARY (NDOF, GK, CLAN, FX, FY, FM, F,BX1,BY1,BM1, BX2, 

& BY2, BM2) 
This subroutine is used to impose the boundary conditions upon 
the global stiffness matrix and force vector. 

INTEGER NDOF, BX1, BY1,BM1,BX2, BY2, BM2,CLAN,i,N,I1,1I2,13,P4 
PARAMETER (P4=99) 
REAL GK (P4,P4),FX, FY, FM, F (P4) 


qaaqaa 


aa 


....invoke the essential boundary conditions...............0005 
IF (BX1 .EQ. 1) THEN 

CALL IMPOSE BC (NDOF,GK,1,F) 
ENDIF 


IF (BY1l .EQ. 1) THEN 
CALL IMPOSE _BC (NDOF,GK,2,F) 
ENDIF 


IF (BM1 .EQ. 1) THEN 
CALL IMPOSE _BC (NDOF,GK,3,F) 
ENDIF 


IF (BX2 .EQ. 1) THEN 

N=NDOF-2 

CALL IMPOSE BC (NDOF,GK,N,F) 
ENDIF 


IF (BY2 .EQ. 1) THEN 

N=NDOF-1 

CALL IMPOSE BC (NDOF,GK,N,F) 
ENDIF 


IF (BM2 .EQ. 1) THEN 
CALL IMPOSE _BC (NDOF,GK,NDOF, F) 
ENDIF 


Cc -...add the concentrated load to the force vector............ 
I1=(CLAN-1) *3+1 
T2=(CLAN~1) *3+2 
I3=(CLAN~1) *3+3 


F (11) *F (I1) +FX 
F(I2) =F (I2)+FY 
F (13) =F (13) +FM 


RETURN 
END 


RRR RK KKKERREKAKEREKRREKKEKEAEREKAEREKREEKEKKKKKKAKKKKKKEKKK 
® * 
SUBROUTINE IMPOSE_BC (NDOF,GK,N,F) 

a Se My ea ar aS ee ee ee 
This subroutine is used to do the redundant leg work of impos- 
ing the boundary conditions. 


aqaagaga 
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100 


...-declare the variables... wc ce cee ee ec wwe cw ewww ce nee 
INTEGER NDOF,N,i,P4 

PARAMETER (P4=99) 

REAL GK (P4,P4),F(P4) 


...-impose the boundary condition on the GK and F arrays...... 
DO 100 i=1,NDOF 
GK(N,i) = 0.0 
CONTINUE 
GK(N,N) = 1.0 
F(N) = 0.0 


RETURN 
END 


HRM IHRE RK REE KKK HK KEE EKER EK KERKKEEREKKEKEKKKKKREAKKEKKEKE 


* 


qaaqaqaana 


* 


SUBROUTINE UPSCALE (NDOF, GK, F, BK, BF) 

This subroutine is used to change the stiffness matrix & force 
vector from single precision to double precision in order to 
solve the linear system of equations in double precision. 

oss A@ClLare the. varlablesi.osoiis is sees FG eG Wd ores 8 wien 8 WHS Siew 6 ae were 
INTEGER NDOF,i, j,P4 

PARAMETER (P4=99) 

REAL GK (P4,P4),F (P4) 

REAL*8 BK(P4,P4),BF(P4) 


-..-generate the doubleprecision compliments of GK and F...... 
DO 110 i=1,NDOF 
DO 100 j=1,NDOF 
BK(i, 3) = GK(i, j) 
CONTINUE 
BF(i) = F(i) 
CONTINUE 


WHR REKKEKKAKKKREKKREKREKREKKKEKKEKRREKKEREKKEKEKEKEKEKKEKKKKEKEK 


* 


aqaaagaagaa 


100 


x 


SUBROUTINE DOWNSCALE (NDOF, BU, U) 

This subroutine is used to do down scale the double precision 
solution of the linear system of equations back to single pre- 
cision. DOT could have problems with double precision numbers! 
.-.-Geclare the variables... ..... cc cece ccc eee eer cece ree ens 
INTEGER NDOF,i,P4 

PARAMETER (P4=99) 

REAL U(P4) 

REAL*8 BU(P4) 


....generate the doubleprecision compliments of GK and F...... 
DO 100 i=1,NDOF 

U(i) = BU(i) 
CONTINUE 


RETURN 
END 


RRR REE KERERKEERREREEKEKKERAEKEKKKEKEKKKKKEKEEKKEKK 
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SUBROUTINE STRESS (X, Y, ALPHA, BETA, U, NEL, ELEN, YOUNG, H, SIGMA_T) 


This subroutine computes the stress at each nodal point. 

see @ARPCTLATACL OMG aia so i5 FS) ce oes aroha 5 Ok ei SES bs, ww ieee Wwe eee Ww wes 

INTEGER NEL,NSNP,NDOF, SWITCH,i,1I1,12,13,14,15,16,17,18,19, 
P1,P2,P4 

PARAMETER (P1=32, P2=33, P4=99) 

REAL ELEN, H,X(P2),Y(P2),ALPHA(P1),BETA(P1), YOUNG, 
K1,K2,CA1,CA2,CB1,CB2,v,vl,v2, 
elxdisp(P2),elydisp(P2),ELEN_f(P2),DISPLEN(P2), 
U(P4),SIGMA_N(P4),SIGMA_B(P4),SIGMA_ T(P4) 


-.--determine the Constants... cece ccc cece cree reece nes 

K1=6.0/ (ELEN**2.0) 

K2=2 .0/ (ELEN) ‘ 
NSNP = (NEL + 1) 

NDOF = NSNP*3 


.-...determine the bending stresses........... cc eee eee eee eee 

DO 100 i=2,NEL 
I1=(i-2) *3+1 
I2=(i-2) *3+2 
I3=(i-2) *3+3 
I4=(i-1)*3+1 
I5=(i-1) *3+2 
I6=(i-1) *3+3 
I7=(i-0) *3+1 
I8=(i-0) *3+2 
I9=(i-0) *3+3 


CBl= COS (BETA (i-1)) 
CAl= COS (ALPHA (i-1)) 
CB2= COS (BETA (i) ) 
CA2= COS (ALPHA (i) ) 


v2 = K1* (U(I4)-U(I1)) *CB1 
+K1* (U(I2) -U(I5) ) *CAl 
+K2* (U(13)+2.0*U (16) ) 


vl = K1* (U(I4)-U(I7) ) *CB2 
+K1* (U(18) -U(I5) ) *CA2 
-K2* (U(I6) *2.0+U (19) ) 


IF (ABS(vl) .GE. ABS(v2)) THEN 
v=vl 
ELSE 
v=v2 t 
ENDIF 
SIGMA_B(i) = YOUNG* (H/2.0) *v 2 
CONTINUE 
v = K1*(U(1)—-U(4) ) *COS (BETA (1) ) 
+K1* (U (5) -U (2) ) *COS (ALPHA (1) ) 
-K2* (U(3) *2.0+U (6) ) 
SIGMA_B(1) = YOUNG* (H/2.0) *v 


v = K1* (U(NDOF-2) -U (NDOF-5) ) *COS (BETA (NEL) ) 
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& +K1* (U (NDOF-4) -U (NDOF-1) ) *COS (ALPHA (NEL) ) 
& +K2* (U (NDOF-3) +2.0*U (NDOF) ) 
Cc 
SIGMA_B(NEL+1) = YOUNG* (H/2.0) *v 
c 
Cc ....determine the normal stresseS...... ee ee ee te eens 
SWITCH = 1 
Cc 
IF (SWITCH .EQ. 1) THEN 
DO 300 i=#2,NEL 
Il = (NEL-2) *3+1 
I2 = (NEL-2) *3+2 
I3 = (NEL-2) *3+3 
I4 = (NEL~-1) *3+1 
P I5 = (NEL-1) *3+2 
I6 = (NEL-1) *3+3 
I7 = (NEL-0O) *3+1 
‘a I8 = (NEL-0) *3+2 
I9 = (NEL-0) *3+3 
Cc 
CA1= COS (ALPHA (NEL-1) ) 
CBl= COS (BETA (NEL-1) ) 
CA2= COS (ALPHA (NEL) ) 
CB2= COS (BETA (NEL) ) 
c 
v2 = (U(I4)-U(I1))*CA1l + (U(I5)-U(I2)) *CB1 
vl = (U(I7)-U(I4))*CA2 + (U(I8)-U(I5) ) *CB2 
Cc 
IF (ABS(vl) .GE. ABS(v2)) THEN 
vevl 
ELSE 
v=v2 
ENDIF 
Cc 
SIGMA_N(i) = (YOUNG/ELEN) *v 
300 CONTINUE 
Cc 
v = (U(4)-U(1))*COS (ALPHA(1)) + (U(5)-U(2)) *COS (BETA(1) ) 
SIGMA_N(1) = (YOUNG/ELEN) *v 
Cc 
v = (U(NDOF-2) -U (NDOF-5) ) *COS (ALPHA(NEL)) + 
& (U (NDOF~-1) -U (NDOF-4) ) *COS (BETA (NEL) ) 
SiGMA_N(NEL+1) = (YOUNG/ELEN) *v 
Cc 
ELSE 
DO 350 i=1,NEL+1 
SIGMA_N(i) = 0.0 
r 350 CONTINUE 
ENDIF 
Cc 
e Cc ....determine the total stresses at each node................. 


DO 400 i=1,NEL+1 
SIGMA_T(i) = ABS(SIGMA_B(i)) + ABS (SIGMA_N(i)) 
400 CONTINUE 


RETURN 


RARER EEEKHRERERAKERAKREKREEREREAREEREEKERKEKEKEKKEKEKKKEKKERKEKEKKEK 


* * 


SUBROUTINE ARCH_OUTPUT 


Cc Re ee ier NGM ES A eye et em TS ee eS ea 
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110 
115 


120 
125 


210 
220 


310 
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This subroutine formats the final results and output of the 
optimization problem and stores it in a file named ARCH_OUT.DAT 
ence tG@ClLAVG” VATTADLOS oe inle cease 6.6 i5p 55 iin we een wi S Mir ae Gide old ee eee tg ee ae 
INCLUDE ‘ARCH_COM.FOR’ 

REAL VOL, VOLUME 





....open output file and write header..............-. cece cence 
OPEN (9, FILE=’ARCH OUT.DAT’, STATUS=’ UNKNO¥N’ ) 


WRITE(9,100) LABEL 
WRITE(9,100) ‘ OPTIMIZATION SOLUTION’ 
WRITE(9,105) ' ----------------~------------------------------ 


FORMAT (/5X, A) 
FORMAT (5X, A) : 
bis ee BOOT L ON A ee at cis fore haved eG a8 ig sos lan ve ale Ye arel osero tela OMe, Biel Shacelev eel ale Biel ee 
WRITE(9,100) ’ A) Problem Parameters: ’ id 
WRITE(9,110) ’ Arch Angle :’, ANGLE, ’ Youngs Modulus:’, YOUNG 
WRITE(9,110) ’ Arch Radius:’, RADIUS, ’ Yield Strength:’,YIELD 
WRITE(9,115) ’ Arch Height:’, H, ‘ No of Elements:’,NEL 
FORMAT (8X,A,F12.3,1T38,A,F12.1) 
FORMAT (8X,A,F12.3,1T38,A,110) 

se 6 SIS OCOCC LON Bis hiner ete aie ene ¥es'o. 55 few siren Olas © Lelie ice: we oh oud lace eee se fare ca seb eee Se 
WRITE(9,100) ’ B) Derived Constants:’ 
WRITE (9,120) ’ No of System Nodal Points...’,NSNP 
WRITE(9,120) ’ No of Degrees of Freedom....’,NDOF 
WRITE(9,125) ‘ Length per Element.......... ', ELEN 
WRITE (9,125) ‘ Phi Angle per Element.......’,PHIANG 
WRITE(9,120) ‘ Number of Iterations........'’,ITERATE 
FORMAT (8X,A, 16) 
FORMAT (8X,A,F12.4) 

2 SBOCTCLON OC aided scale Ge Gidie Sib SS alo eos i bw eat aie ae Wcaierie: eke el whale So ecer arate ce 
WRITE (9,100) ’ C) Structure Loading:’ 
WRITE (9,125) EX ots. cece eveeseewas 04 ee release ’,FX 
WRITE (9 7:125) EX oo iio: oselete sre ovale G.0. ies ehe eee ‘,FY 
WRITE (9,125) ‘EM... cc cece twee ec cence ‘,FM 
WRITE (9;,.225) RA bo iceks 066.6 505 s.e eee So we ee eee ‘,FA 

See SSCELOM Ds oe snes ei ai iee Sheer obi Sse isi er a iéiere: 6 eoeie veces id lore le, owes order Sie Seiiaire 
WRITE(9,100) ’ D) Elemental Dimensions and Stress Distribution:’ 
WRITE (9,210) ‘Element’, ’Height’,’Base’,’Length’,’ Volume’ 
FORMAT (8X,A,T19,A,1T34,A,T50,A, T62,A) % 
FORMAT (8X, 14,T17,F10.5,1T32,F10.5,7T48,F8.5,T60,F8.5) 
VOLUME = 0.0 

8 


DO 300 i=1,NEL 
VOL = H*ELEN*BASE (i) 
WRITE (9,220) i,H, BASE (i) ,ELEN, VOL 
VOLUME = VOLUME + VOL 

CONTINUE 


Silas CBECEL ON El ese ccd we Oe ee ie aS wha en aad taal te crear ee eet are alse 
WRITE(9,100) ’ E) Objective Function:’ 

WRITE (9,310) ’ Total structure Volume:’ , VOLUME 

FORMAT (8X,A,F12.6) 
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320 
330 


400 
410 
420 
430 


WRITE (9,330) ‘’Node’,’Stress’ 
DO 320 i=1,NSNP 

WRITE (9,*) i,SIGMA_T(i) 
CONTINUE 
FORMAT (8X,A,T19,A) 


Sie sSSCCE LON: ME ies ss aiooieve ee ees Sanh S Ge wee whee Coe dren duleew ee o-eiy ere eee RS 


WRITE(9,100) ’ F) Boundary Conditions:’ 


WRITE(9,410) ‘Node’,’X-Displ’,’Y-Displ’,’ Slope’ 


WRITE (9,430) 1,BX1,BY1,BM1 
WRITE (9,430) NEL+1,BX2, BY2, BM2 


Sb ge BO CE LON MG a or cidovend Sie Scel teh helo sens S Sek eS wg Owen a btw e Sud-dis ele Biel. 4 eke 


WRITE(9,100) ’ G) Solution Vector:’ 


WRITE(9,410) ‘Node’,’X-Displ’,’Y-Displ’,’Slope’ 


DO 400 i=1,NSNP 
I1*(i-1) *3+1 
I2=(i-1) *3+2 
I3=(i-1) *3+3 
WRITE (9,420) i,U(I1),U(12),U(13) 
CONTINUE 
FORMAT (T9,A,T17,A,T31,A,T46,A) 
FORMAT (7X, 15, 3E14.6) 
FORMAT (7X, 1I5,7T20,14,1T34,14,1T48,14) 


RETURN 
END 
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Cc ARCH_COMMON 
c 
Cc oe. et SAS@ELN ECL ONS ores 5 sec were) crea ae oe Sere ale SINS Seas Shae O000S gee Oe a Sees 
Cc Blo. sac The maximum number of elements 
Cc B2 esses The maximum number of global nodal points 
Cc |< err The maximum number of design constraints 
c PAs feeds The maximum number of degrees of freedom 
Cc 
Cc oeeeGeclare the variables. so. ccc ceed ce cats Saeae eee ee seewes 
INTEGER NEL, NCON, NSNP, NDOF, METHOD, MINMAX, INFO, IPRINT, IWK (400), 
& NRWK, NRIWK, IPRM(20) , COUNT, OPTDCS, ITERATE, PRCSN, CLAN, 
é BX1, BY1,BM1,BX2,BY2,BM2,P1,P2,P3,P4 
Cc 
PARAMETER (P1=32, P2=33, P3=96, P4=99) 
Cc 
REAL ANGLE, RADIUS, ELEN, H,X(P2),Y(P2),ALPHA(P1),BETA(P1), 
& YOUNG, YIELD, WK (27000) , RPRM(20) ,OBJ,G(P3), 
é DV1BG, DV1LO, DV1UP, BASE (P1), BASEL (P1),BASEU(P1), 
& FA, FX, FY, FM, U(P4) , SIGMA_T (P4) ? 
Cc 
c aot MARS AN: COMMON 2:6 6255550. co: 6:10 lee Sal Sie le See la Pek sega eile eretele e w-a. 9S) 0 4S ose S 
COMMON NEL,NCON,NSNP,NDOF, METHOD, MINMAX, INFO, IPRINT, IWK, 
& NRWK, NRIWK, IPRM, COUNT, OPTDCS, ITERATE, PRCSN, CLAN, 
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